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We introduce a new kludge scheme to model the dynamics of generic extreme mass-ratio inspirals 
(stellar compact objects spiraling into a spinning supermassive black hole) and to produce the gravi- 
tational waveforms that describe the gravitational-wave emission of these systems. This scheme com- 
bines tools from different techniques in General Relativity: It uses a a multipolar, post-Minkowskian 
expansion for the far-zone metric perturbation (which provides the gravitational waveforms, here 
taken up to mass hexadecapole and current octopole order) and for the local prescription of the 
self- force (since we are lacking a general prescription for it); a post-Newtonian expansion for the 
computation of the multipole moments in terms of the trajectories; and a BH perturbation theory 
expansion when treating the trajectories as a sequence of self-adjusting Kerr geodesies. The orbital 
evolution is thus equivalent to solving the geodesic equations with time- dependent orbital elements, 
as dictated by the multipolar post-Minkowskian radiation-reaction prescription. To complete the 
scheme, both the orbital evolution and wave generation require to map the Boyer-Lindquist coordi- 
nates of the orbits to the harmonic coordinates in which the different multipolar post-Minkowskian 
quantities have been derived, a mapping that we provide explicitly in this paper. This new kludge 
scheme is thus a combination of approximations that can be used to model generic inspirals of sys- 
tems with extreme mass ratios to systems with more moderate mass ratios, and hence can provide 
valuable information for future space-based gravitational-wave observatories like the Laser Interfer- 
ometer Space Antenna and even for advanced ground detectors. Finally, due to the local character 
in time of our multipolar post-Minkowskian self-force, this scheme can be used to perform studies 
of the possible appearance of transient resonances in generic inspirals. 



PACS numbers: 04.25.Nx,04.25.-g,04.30.Db,04.30.-w 

I. INTRODUCTION 

Gravitational waves (GWs) hold the promise to pro- 
vide detailed information about astrophysical bodies that 
are obscure in the electromagnetic spectrum, such as bi- 
nary black hole (BH) systems. Moreover, such waves 
will allow for the first studies of the nature of the grav- 
itational interaction and of the validity of General Rel- 
ativity (GR) in the strongest regimes [H [2]. An accu- 
rate modeling of such GWs is essential for the extraction 
and characterization of weak signals buried in detector 
noise. This is because waveform templates act as an op- 
timal linear filter that maximizes the signal-to-noise ratio 
(SNR) in the presence of stochastic noise. The absence of 
such templates for certain GW sources renders subopti- 
mal any GW search strategy. Therefore, the construction 
and modeling of GWs to construct accurate templates for 
data analysis is of paramount importance in the blossom- 
ing of GW astrophysics. 

One of the staple GW sources of the planned space- 
based observatory Laser Interferometer Space Antenna 
(LISA) [31 |3] that require accurate templates for their 
detection and analysis are extreme-mass-ratio inspirals 
(EMRIs) [5]. These events consist of a small compact 
object (SCO), such as a stellar-mass BH or neutron star 



(NS), spiraling in a generic orbit into a spinning, (su- 
permassive black hole (MBH), whose evolution is GW 
dominated. In such inspirals, the SCO spends up to mil- 
lions of cycles in close orbits around the MBH, possibly 
with large pericenter velocities and eccentricities, sam- 
pling the strong gravitational field of the MBH. 

Many astrophysical scenarios predict the existence of 
such EMRIs. One such scenario postulates that the SCO 
exchanges energy and angular momentum with other 
stars in a stellar core/cusp near a MBH at a galactic cen- 
ter, via two-body relaxation and dynamical friction [5]. If 
so, the SCO can be swung sufficiently close to the MBH 
to be gravitationally captured [5, 6 , at which point it 
would slowly inspiral until being swallowed by the MBH. 
Of course, such a scenario is complicated by mass seg- 
regation [7], triaxial density profiles [8], resonant relax- 
ation [5], etc. Other channels of EMRI formation include 
binary tidal separation [TU] (where a binary is disrupted 
with one component captured by the MBH) and mas- 
sive star capture or production in accretion discs [llj . 
where the stellar-mass BH is directly formed in the ac- 
cretion disk of the MBH. It turns out that the EMRI 
GWs produced by EMRIs in each astrophysical mech- 
anism have distinct characteristic that may be use to 
distinguish them from EMRI observations [T2] . 
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The inspiral is by far the dominant GW phase for 
EMRI data analysis purposes. This can be understood 
rather easily (see, e.g. [T3J Q3]). The number of cycles 
accumulated in the inspiral scales with the inverse of the 
mass ratio ro*/M,, while in the plunge and merger it 
scales with the MBH mass M, only. Since the mass ratio 
for EMRIs is in the range 0(l(r 4 )-0(lCr 6 ), the charac- 
teristic duration and cycle accumulation during the in- 
spiral phase is several orders of magnitude larger than 
during plunge and merger. In turn, the SNR is in the 
range O(10)-O(10 2 ) for EMRIs at realistic distances and 
it scales linearly with the total number of cycles. Thus, 
the contribution of the plunge and merger to the SNR 
is reduced by a factor of 0(m*/M,) relative to the in- 
spiral contribution. In addition, there is no detectable 
ringdown for EMRIs, as the SCO barely perturbs the 
background geometry as it crosses the MBH's event hori- 
zon. Since the SCO is not disrupted in EMRI plunges, 
the SCOs internal structure is erased or effaced almost 
completely, without affecting the inspiral signal. 

A positive consequence of the large number of EMRI 
GW cycles is that these waves carry a detailed map of 
the MBH geometry, so that we expect to determine the 
EMRI physical parameters with high precision [TS] . This 
information will be very useful, in particular to test the 
spacetime geometry of MBHs jTBl HZ] and even alter- 
native theories of gravity (see, e.g. [T7H2D]). Moreover, 
given that the expected even rate is in the range 10 — 10 3 
EMRIs/yr [3,(21], EMRI observations will allow us to un- 
derstand better the stellar dynamics near galactic nuclei, 
populations of stellar BHs, etc. (for a review see [5]), and 
it also possible that they will tell us about cosmology [55] . 

Although one is left to model the inspiral phase of EM- 
RIs, this task remains a gargantuan endeavor for several 
reasons. First, the accuracy requirements for EMRI tem- 
plates are much more stringent than for comparable-mass 
binaries. For detection and parameter estimation one 
usually demands an absolute accuracy of better than 1 
and SNR -1 radians in the GW phase respectively over 
the entire time of observation. In a 1 yr inspiral, a 
typical EMRI can have 10 5 cycles in-band, which then 
translates into a relative radian accuracy of 0(1O -6 ) and 
0(1O -8 ) for detection and parameter estimation respec- 
tively. It is important to note that these precision re- 
quirements are just simple estimates that do not take 
into account data analysis strategies that may relax them 
(see, e.g. [5T]). In contrast, the above relative measure 
becomes of O(10 _2 /SNR) for ground-based data anal- 
ysis of comparable-mass plunge-merger-ringdowns with 
current numerical relativity simulations, because these 
accumulate only O(10) GW cycles in these phases. In be- 
tween EMRIs and comparable-mass inspirals, there are 
Intermediate-Mass-Ratio Inspirals (IMRIs) , which can be 
potential sources for both space-based detectors (where 
an Intermediate-Mass BH (IMBH) inspirals into a MBH) 
and advanced ground-based detectors (where a SCO in- 
spirals into an IMBH). 

The extreme mass ratios involved in the problem also 



lead to the appearance of two different spatial and time 
scales. The two different spatial scales are represented 
by the very different sizes of the MBH and the SCO, 
m*/M. <C 1, whereas the different time scales are the 
orbital one and the one associated with the radiation- 
reaction effects, T 0lbital /T Rn ~ m+fM, -C 1. An illustra- 
tion of how this complicates the EMRI problem is the 
recent work of Lousto and Zlochower [53], who evolved 
the first 1 : 100 mass-ratio binary over the last two orbits 
before merger and plunge; this simulation took approxi- 
mately 100 days of computational time using full numer- 
ical relativity. This means that with present numerical 
relativity techniques, full numerical simulations are out 
of the question for EMRI modeling. 

Another key reason for the difficulty of EMRI mod- 
eling is their intrinsic strong-relativistic nature. In the 
interesting part of the EMRI dynamics, the SCO is mov- 
ing around the strong-field region of the MBH, acquiring 
large pericenter velocities and even sampling regions in- 
side the ergosphere, leading to large relativistic T-factors 
over tens of thousands of GW cycles. Approximate tech- 
niques employed in the comparable-mass regime, such as 
low- velocity expansions in the post-Newtonian (PN) ap- 
proximation, are then ill-suited for EMRIs. So neither 
numerical relativity nor PN theory arc, for very different 
reasons, suitable schemes to model EMRIs. 

A better-suited framework that exploits the extreme 
mass ratios involved is BH perturbation theory, where 
one treats the SCO as a small perturbation of the MBH 
background geometry. In this context, the inspiral can be 
described as the action of a self-force. This local vector 
force is made out of the regularized metric perturbations 
generated by the SCO, after eliminating divergences due 
to the particle description of the SCO. The SCO's motion 
is then governed by the MiSaTaQuWa equation of mo- 
tion equation, derived in J24[ [25] (see also Sec. [TTJ) . The 
MiSaTaQuWa equation is considered the foundation of a 
self-consistent scheme to describe EMRIs in an adiabatic 
way by coupling it to the partial differential equations 
that describe the perturbations produced by the SCO. 
For recent discussions on these issues see [5"6"H55] and for 
general reviews see [5M3T] . 

At present, the gravitational self-force has been com- 
puted for the case of a nonrotating MBH using time- 
domain techniques [35] (see also [35H35] for the study of 
the physical consequences of the self-force) and progress 
is being made towards calculations for the more astro- 
physically relevant case of a spinning MBH [36]. In the 
meantime, a number of new techniques in the frequency 
and time domains are being developed to produce ac- 
curate and efficient self- force calculations [37]. In any 
case, given the amount of cycles required for EMRI GWs 
and the present complexity of self-force calculations, we 
cannot expect to generate complete waveform template 
banks by means of full self-force calculations. Instead, 
the goal of these studies should be to understand all the 
details of the structure of the self-force so that we can 
formulate efficient and precise algorithms to create the 
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waveforms needed for LISA data analysis, perhaps com- 
plementing some of the existent approximating schemes 
that we review below. 



A. Existent EMRI Waveform Models 

In parallel to the efforts to make progress in the 
self-force program, there have also been some efforts 
to build certain less-reliable approximation schemes to 
model EMRIs. These are very useful, for instance, for 
parameter estimation studies. We compare and contrast 
these in Table [TJ ordered by level of complexity from sim- 
plest (top) to most complex (bottom), where we also in- 
clude the new Kludge scheme at the bottom for compar- 
ison. 

The simplest model is that of Peters and Mathews [38] , 
in which the SCO is assumed to be moving on Ke- 
plerian ellipses. The orbital elements of this ellipse 
evolve according to leading-order (Newtonian), dissi- 
pative radiation-reaction, i.e. the quadrupole formula 
for the loss of energy and angular momentum. Wave- 
forms are then computed also to leading order via the 
quadrupole formula |44j . 

A better model was introduced by Barack and Cut- 
ler |15j , the so-called Analytical Kludge waveform model. 
This model is based on the Peters & Mathews 38 model, 
but it is enhanced via different PN formulas in order to 
account for all the relativistic effects, both dissipative 
and conservative, present in a generic EMRI event. It 
has the advantage that the orbital evolution is decoupled 
from the evolution of the additional relativistic effects 
(which evolve effectively in the radiation-reaction time 
scale), and hence EMRI waveforms, also prescribed via 
the quadrupole formula [44j . can be computed very fast. 
For this reason, it has become the method of choice for 
many LISA parameter estimation and data analysis stud- 
ies (see, e.g. [4"5]). 

A more sophisticated approach was proposed by 
Babak, et. al. [39] and is sometimes referred to as the 
Numerical Kludge waveform model. In this setup, the 
orbital motion is given by geodesies around a Kerr BH 
and the radiative effects are prescribed via PN evolution 
equations for orbital elements (from 2PN expressions for 
the fluxes of energy and angular momentum) calibrated 
to more accurate Teukolsky fluxes with 45 fitting param- 
eters |46j . The waveforms are then modeled again via a 
multipolar expansion [47], but this time taken to next- 
to- leading-order (quadrupole plus octopole). 

Recently, a new hybrid scheme has been proposed by 
Yunes, et. al. [HI HOI IUJ based on effective-one-body 
(EOB) techniques [48] . In this approach, the SCO- 
MBH, two-body system is mapped to an effective one- 
body system: a Kerr BH perturbed by a small effective 
object. The orbital motion is obtained by solving the 
Hamilton equations for the Hamiltonian of the effective 
system. When neglecting conservative self-force correc- 
tions, this reduces to solving the geodesic equations in the 



Kerr background. The radiation-reaction comes from the 
short- wavelength approximation of Isaacson's |49j , where 
the waveform is constructed from an orbit-averaged, but 
resummed PN expression |50j . This radiation-reaction 
force is then enhanced through the addition of very high 
PN order point-particle results [UJ [52] and expressions 
that account for the flux of energy and angular momen- 
tum into the MBH's horizon [5U[S3]. Although shown to 
be accurate for equatorial, circular orbits ED] relative to 
Teukolsky waveforms, the EOB scheme has not yet been 
tested for eccentric or inclined orbits. 

Another approach is the Teukolsky-b&sed scheme of 
Hughes and others [121 331 |54"H6"3"] have developed. In 
this model, one prescribes the inspiral as a sequence of 
slowly-changing geodesies. The mapping between them is 
given by the orbit-averaged evolution of orbital elements, 
which in turn is obtained by a balance law relating av- 
eraged fluxes at the boundaries of spacetime and at the 
location of the SCO. These averaged fluxes at a given 
geodesic or point in orbital phase space are computed 
by solving the Teukolsky equation at that point. There- 
fore, the construction of any single waveform requires the 
mapping of the entire orbital phase space, which in turn 
is computationally prohibitive for truly generic EMRIs. 
Moreover, this scheme has numerical (and conceptual) 
difficulties when modeling EMRIs in regimes of space- 
time where the evolution deviates from adiabaticity, such 
as in or close to plunge, around rapidly spinning MBHs 
(a/M. > 0.9), or highly eccentric inspirals. 



B. The New Kludge Scheme 

In this paper, we devise a new kludge approxima- 
tion scheme (relative to numerical kludge) that combines 
seemingly disparate ingredients from BH perturbation 
theory and the multipolar post-Minkowskian formalism 
of Blanchet and Damour [64] . We shall here combine two 
different approximation schemes, BH perturbation the- 
ory (which assumes only that the mass-ratio is small) and 
the multipolar post-Minkowskian formalism (which as- 
sumes the gravitational field strength is small), together 
with other ingredients related to the choice of coordinate 
system and waveform construction method. 

In the new kludge scheme, the orbital motion is 
prescribed as a spacetime trajectory that is piecewise 
geodesic (with respect to the MBH geometry) and such 
that the different geodesic intervals are connected via the 
SCO's local, self-acceleration (due to its own gravity in 
the presence of the MBH). In this sense, each geodesic 
interval can be chosen arbitrarily small. This is in con- 
trast to the Teukolsky approach (see, e.g. [42]) where 
the mapping between sequences is given by the averaged 
(over several orbits) GW energy-momentum fluxes and 
balance laws. This fact, that the mapping is given purely 
in terms of quantities local to the SCO's worldline and 
not by nonlocal balance laws, is a distinctive feature of 
the new kludge scheme. The implementation of this idea 
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Scheme Name 
Peters & Mathews [38f 
Analytic Kludge fTS] 
Numerical Kludge [39] 

EOB [I4ll40l l4l] 
Teukolsky-based [4*2"1H5] 
This work 



Orbital Motion 



Radiation- Reaction 



Waveform Generation 



Keplerian Ellipses 
Keplerian Ellipses 
Kerr Geodesies 
Kerr Geodesies 
Kerr Geodesies 



Newtonian Order 
Low-Order PN 
Calibrated Low-Order PN 
Calibrated, Resummed 5.5PN 
Averaged Teukolsky 



Kerr Geodesies Local Multipolar Post-Minkowskian 



Multipolar decomposition (I = 2) 
Multipolar decomposition (I — 2) 
Multipolar decomposition (I < 3) 
Resummed 5.5PN (I < 8) 
Adiabatic Teukolsky (I < 60) 
Multipolar decomposition (I < 4) 



TABLE I. Comparison of existent modeling schemes. 



uses the evolution of geodesies with varying orbital el- 
ements; the energy, angular momentum in the spin di- 
rection, and Carter constants are then functions of time 
governed by the self-force. 

Such an evolution scheme is a direct implementation 
of the osculating orbits method, proposed by Pound and 
Poisson [65] and Pound [27] . This method professes that 
at each point of the SCO's nongeodesic worldline there is 
a unique geodesic that lies tangent to it. Therefore, the 
worldline is simply an interpolation between these tan- 
gent geodesies. Such a scheme hinges on a fundamental 
assumption of EMRI modeling: the adiabatic approxima- 
tion, which assumes the deviation vector between adja- 
cent tangent geodesies is small, i.e. the radiation-reaction 
time scale is much longer than all other timescales, par- 
ticularly the orbital one. As Gralla and Wald [53] ex- 
plained, the adiabatic approximation only holds quasilo- 
cally around some small neighborhood of proper time at 
each point of the SCO's worldline. This problem can be 
circumvented, however, if at each point on the worldline 
the deviation vector is recomputed, which is exactly the 
basis of the method of osculating orbits [27l [65] . 

In the new kludge scheme, the evolution of orbital el- 
ements is prescribed by the SCO's self-acceleration, but 
this quantity is not known exactly (numerically or other- 
wise) for generic orbits around a spinning MBH. In view 
of this, we model the self-force via a radiative approx- 
imation, i.e., through the time-asymmetric part of the 
radiation field, given by the "half-retarded minus half- 
advanced" Green function |66j . This can be implemented 
with a high-order multipolar, post-Minkowskian expan- 
sion, which assumes gravitational radiation is small rel- 
ative to the gravitational field of the background, i.e. an 
expansion in powers of Newton's gravitational constant 
G. 

The new kludge, local self-force prescription is com- 
pleted once we say how the multipole moments depend on 
the orbital trajectories. This can be achieved by asymp- 
totically matching a PN and a post-Minkowskian solu- 
tion [6TH72) . In this paper, we use only the leading-order 
expressions for these multipoles, although in the future it 
would be trivial to including higher PN corrections. As 
such, we are not consistently keeping a given PN order, 
but instead using leading-order expressions for the mul- 
tipole moments, while keeping several such moments in 
the expansions. 

The use of the radiative, multipolar post-Minkowskian 
expansion is only because of the lack of a more precise 



self-force. Clearly, this radiative approximation neglects 
the conservative part of the self-force, which could be 
important in GW modeling [65]. Once the full self- force 
becomes available, however, one could easily employ it 
instead of its post-Minkowskian expansion. Our set-up 
is general and easily adaptable to other, more precise 
expressions for the self-force. 

Once the orbital evolution has been prescribed, one 
can construct the GWs again through multipolar, post- 
Minkowskian expressions in terms of a sum over mul- 
tipole moments. Since the mapping between Boyer- 
Lindquist and harmonic coordinates is known, there are 
no coordinate issues to relate the trajectories obtained 
from the orbital evolution to the trajectories that enter 
the definition of the multipole moments. We here employ 
an expansion to second-order in the multipole moments, 
including both the mass hexadecapole and the current oc- 
topole, thus keeping contributions one order higher than 
traditional kludge waveforms. 

Let us emphasize that the new kludge waveforms are 
only approximate gravitational-wave solutions, meant to 
be used for qualitative descoping studies and investi- 
gations of resonant behavior. That is, the waveforms 
constructed could be used to determine the accuracy to 
which parameters could be extracted given a detection 
with new space-based gravitational- wave observatories as 
a function of the detector. Moreover, one could also study 
how the resonances found by Flanagan and Hinderer [75] 
affect such parameter estimation and detection. Having 
said that, improvements would have to be implemented 
before our kludge waveforms can be used as realistic tem- 
plates in data analysis. As we will see, the new kludge 
scheme is amenable to such improvements, which will be 
the focus of future work. 



C. Comparison to Standard Kludge Waveforms 

The main advantage of the new kludge scheme is 
the fact that the radiation-reaction effects are described 
in terms of a local self-force. This means that the 
inspiral description does not need to average certain 
gravitational- wave fluxes over a number of periods/orbits 
like is traditionally done in kludge implementations. In- 
stead, the self-force is prescribed through a multipolar, 
post-Minkowskian expansion (e.g. the quantity A Rn men- 
tioned above) that contains time-derivatives of the sys- 
tem multipole moments in a nonaveraged form. 
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When implementing such a nonaveraged and local 
scheme, it is critical to use an exact mapping between 
Boyer-Lindquist coordinates, used in the integration of 
the geodesic equations of motion, and harmonic coor- 
dinates, employed in the calculation of the multipolar, 
post-Minkowskian self-force. This eliminates gauge is- 
sues that plague kludge waveforms due to the neglect of 
such a mapping (i.e. kludge schemes simply use Boyer- 
Lindquist like Cartesian coordinates in the multipolar de- 
composition of the GWs). 

The use of a local self-force provides the freedom to 
choose how often to apply radiation-reaction effects in the 
numerical implementation of the dynamics (trajectory 
and waveform construction) of the new kludge scheme. 
The two extremes are: (i) We can apply the self-force 
at every single time step, which corresponds to the case 
of a continuous local self-force, or (ii) we can store the 
information about the self-force and apply it after a cer- 
tain period of time, mimicking the averaging procedure 
of other schemes. 

These two extreme ways of using the new kludge 
scheme are in correspondence with two very relevant po- 
tential applications. Whereas the type of application 
(ii) can be used to try to generate efficiently EMRI 
gravitational-wave templates for parameter estimation 
and data analysis development purposes, the type of ap- 
plication (i) seems to be very well fitted for studying local 
phenomena in the dynamics of EMRIs. In this sense, an 
important application of our scheme would be to study 
the transient resonances that Flanagan and Hinderer [73] 
have recently reported in generic EMRI orbits (eccentric 
and inclined) when the fundamental orbital frequencies 
become commensurate. The new kludge scheme can be 
in principle used to study such behavior and shed some 
light on the relevance that it may have for future EMRI 
detection with LISA-like detectors. 

Let us clearly state, however, that the new kludges 
presented here are not intended to model comparable- 
mass systems, for which the effective-one-body frame- 
work has proven the most successful |48j . Instead, we 
here focus on EMRIs only and we are interested in study- 
ing local phenomena that no other kludge scheme can 
currently handle. Perhaps in the future, one could also 
use effective-one-body tools for such studies, but this 
would require a tested framework capable of handling 
completely generic orbits. Work along this lines is cur- 
rently underway [T4l |40| ITT] . 



D. Executive Summary of Main Results 

In this paper we present the new kludge scheme, intro- 
ducing one by one each of its ingredients, from the form 
of the equations of motion for the inspiral to the GW 
construction, including the details of the approximations 
that we use to construct the local self-force that drives 
the inspiral. From a technical point of view, one of the 
main challenges of the numerical implementation of our 



scheme is the computation of high-order time derivatives 
(of the mass and current multipole moments), which are 
crucial for the estimation of the radiation-reaction effects 
(the multipolar post-Minkowskian self-force involves up 
to eight-order time derivatives of the trajectory in har- 
monic coordinates) and the GW construction (since we 
are using up to the mass hexadecapole and current oc- 
topole multipoles in the calculation, we require up to 
fourth-order time derivates of the trajectory in harmonic 
coordinates). The computation of these time derivatives 
is very challenging, forcing us to implement numerical 
techniques adapted to the properties of EMRIs dynam- 
ics. The key point is to use the fact that geodesic orbits 
have, in the generic case of eccentric and inclined tra- 
jectories, three fundamental frequencies. This then al- 
lows us to fit any quantity that needs differentiating to 
a multiple Fourier series, using a standard least-squares 
technique. Numerical derivatives of such quantities can 
then be obtained simply by analytically differentiating 
the Fourier expansion. Numerical experimentation has 
shown that this technique works remarkably well, even 
for the highest-order time derivatives that the new kludge 
requires. 

To illustrate the capabilities of our scheme we show in 
this paper results from evolutions for different types of 
orbits: circular equatorial, eccentric equatorial, circular 
inclined, and the generic eccentric inclined orbits. Using 
these evolutions we study different aspects of the new 
kludge scheme. First, we consider the impact of harmonic 
coordinates in the trajectories and waveform observables 
in comparison with using other coordinate systems. We 
find that not properly accounting for this transformation 
can lead to huge errors in the amplitude and phase of the 
waveform, eg. up to errors of order a factor of 2 in the 
total accumulated cycles after a 1 yr evolution. Second, 
we study the impact of the different radiation-reaction 
potentials in the resulting waveforms. Although these 
corrections have a smaller impact than the proper use of 
coordinates, they are still large for strong-field EMRIs. 
For example, including higher-order terms to the Burke- 
Thornc potential leads to corrections of order 10 4 radians 
in a two month evolution. Third, we investigate the use 
of a quadrupole waveform prescription versus a more ac- 
curate hexadecapole-octopole prescription. We find no 
change in the resulting waveform phases, but an ampli- 
tude correction of less than 5%. Fourth, we consider the 
time-evolution of different orbital parameters when we 
apply radiation-reaction effects locally in time through 
our multipolar, post-Minkowskian self-force. Although 
there are some orbital parameters whose time evolution 
is not affected, other parameters like the inclination an- 
gle can acquire small oscillations with period equal to the 
orbital period. Such small oscillations are not captured 
if the self-force is orbit-averaged. 

Finally, we perform some tests and comparisons with 
results in the literature to validate the new kludge numer- 
ical implementation. In particular we test the prediction 
that the inclination angle remains almost constant during 
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the evolution, a test that our scheme successfully passes. 



E. Notation and Organization of the Paper 

Throughout this paper we use the metric signature 
(—,+,+,+) and geometric units in which G = c = 1. 
The MBH geometry, whose metric we denote by g*a, is 
determined by its mass M. and (magnitude of the) spin 
angular momentum S t = M.a, with dimensionless, Kerr 
spin parameter a/M, (-1 < a/M, < 1). The SCO is 
parameterized only in terms of its mass to* since we ne- 
glect its spin and other internal properties. The binary 
system's parameters are the mass ratio q = m^/M, and 
the total mass M Tot = m* + M,. The reduced mass of 
the system is therefore \i = m*M,/M Tot , while the sym- 
metric mass ratio is r] = fi/M Tot = q/(l + q) 2 . The mass 
difference is denoted by 6m = M, — to*. 

The SCO orbit can be parameterized in terms of the 
constants of geodesic motion I A — (E,L Z ,Q/C), which 
stand for the SCO's energy normalized with respect to 
to*, the z-component of the angular momentum normal- 
ized to to*, and the Carter constant also normalized to 
to* (C and Q stand for two definitions of the Carter con- 
stant that we use in this paper). Alternatively, we will 
also parameterize the SCO orbit in terms of the orbital el- 
ements O a = (e,p, t/0 inc ), which are also constants of the 
geodesic motion, where e is the orbit eccentricity, p is the 
semilatus rectum, and i and 6* inc are two measures of the 
orbit inclination. We present the mappings between X A 
and O a in Appendix |e| The SCO spacetime trajectory 
is denoted via z^(t), where r is proper time and thus, its 
four-velocity is the unit timelike vector u M = dz^/dr. 

Post-Newtonian orders always refer to a relative or- 
dering scheme (instead of an absolute one), such that 
the N-th PN order term refers to one of the form A = 

^Newtonian [! + ■ • ■ + O (v 2N / C 2N )} , where ^Newtonian IS the 

leading order contribution. We shall commonly drop the 
factor of l/c when referring to PN expansions. 

Greek letters in index lists are used to denote indices 
on the 4-dimcnsional spacetime, while Latin letters in 
the middle of the alphabet k, . . . denote spatial in- 
dices only. Covariant differentiation is denoted using the 
symbol V„, while partial derivatives with respect to the 
coordinate x M are denoted as d p B v or B . We denote 
symmetrization and antisymmetrization with parenthesis 
and square brackets around the indices respectively, such 

as \u) = IV + A vM 2 and a \m = IV - A vM 2 - 

We use two main sets of coordinate systems: Boyer- 
Lindquist coordinates Xb L and harmonic coordinates 



Xh- Other systems of coordinates that we also consider 
are asymptotic-Cartesian mass-centered (ACMC) coor- 
dinates, Xacmcj an d approximate harmonic coordinates, 
a; Retarded time is denoted in harmonic coordinates 



via i_ 



where r H 



zl) 1 ' 2 . When 



we refer to the Kerr metric, we sometimes use the label 
K, e.g. g*a and g"^. In some situations, it is crucial to 
specify in which coordinate system the metric has to be 



written in a certain equation, like in a coordinate trans- 
formation. In those situations, we also incorporate in 
the metric (or related objects) a label associated with 
the coordinate system, e.g. g a 'p or g„^ L - As for an- 
gular coordinates, we shall find it convenient to some- 
times perform multipolar decompositions as in |47j , with 
spin-weighted spherical harmonics __ 2 Y im and symmet- 
ric and trace-free spherical harmonic tensors y^, where 
L = (ix, i2, . . . , i n ) is a multi-index (see |47j for details 
on the multi- index notation). 

The organization of this paper is as follows. Section [TT] 
introduces the self-force approach to EMPJ modeling and 
some basic details of the method of osculating orbits. 
Section [ill] describes in detail the new kludge scheme, 



including the MBH geometry and the properties of the 
geodesic orbits, the application to them of the method of 
osculating orbits, the multipolar post-Minkowskian ap- 
proximation to the self-force we use and the radiation- 
reaction potentials from which it can be obtained, the 
mapping between Boyer-Lindquist and harmonic coordi- 
nates, and the multipolar post-Minkowskian waveform 
generation formalism that we use in our scheme. Sec- 
tion |IV| explains the numerical implementation of the new 
kludge approach and the main numerical techniques that 
we use. This includes the algorithms for the integration 
of the ODEs governing the local geodesic motion and 
the accurate estimation of time derivatives of several or- 
ders of the multipole moments. Section [V] summarizes 
the different ways in which our scheme can be used and 
presents several numerical results that illustrate its main 
features, using from circular-equatorial orbits to eccentric 
inclined orbits. Section IVTl concludes and discusses sev- 
eral avenues for future research applying the new kludge 
scheme. 

We have attempted to present the main ingredients 
of our approach in the main body of the paper, rel- 
egating some details to the Appendices. Appendix [A] 
gives explicit expressions of the different pieces of the 
multipolar post-Minkowskian self-force in terms of the 
radiation-reaction and local potentials. It also gives for- 
mulas to simplify the computation of the different spatial 
derivatives of the Kerr local potentials. Appendix [B] pro- 
vides complementary formulas related to the mapping 
between Boyer-Lindquist and harmonic coordinates. In 
particular, we give expressions for the components of the 
Jacobian and Hessian, and for the components of the 
covariant and contravariant Kerr metric tensor in har- 
monic coordinates. Appendix [C] constructs a system of 
asymptotically, mass-centered coordinates and a system 
of approximate harmonic coordinates that we compare 
with the exact harmonic coordinates of Sec. IIID| Ap- 



pendix [D] performs a far-field expansion of the Kerr met- 
ric coefficients in the approximate harmonic coordinates 
of Appendix [C] Appendix [E] describes in detail how to 
implement the one-to-one mapping between the orbital 
elements O a = (e,p, i/0 iac ) and the constants of motion 
T A = (E,L Z ,C/Q). Appendix [F] summarizes the main 
formulas for the computation of the fundamental frequen- 
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cies and periods with respect to the Boyer-Lindquist co- 
ordinate time. Finally, Appendix [G| provides expressions 
for the coefficients that determine the evolution of the ra- 
dius and Carter constants of an inspiral through circular 
nonequatorial geodesies. 

II. THE SELF-FORCE APPROACH TO EMRIS 

In this section we review some of the basic and well- 
established concepts related to the self-force as they are 
related to the new kludge approach (see the review pa- 
pers [301 [3H [74] for details). 

The foundations for the first-order perturbative de- 
scription of EMRIs were laid down in the papers by Mino, 
Tanaka, and Sasaki [24] and Quinn and Wald [25 . The 
main result of these papers was the equation of motion 
of a massive pointlike object (describing the SCO) in the 
geometry of a MBH. The stress-energy tensor of the SCO 
is then given by 

-^L^-Z^T)]^-^-, (1) 

where g denotes the determinant of g Q/3 . Then, the SCO 
generates metric perturbations, h a ^, around the MBH 
background geometry that in the Lorenz gauge, 

v,r = o, h^ = h^-\g a ^h, h = g ^h^, (2) 

satisfy the linearized Einstein equations 

Dh a0 + 2R a fl ^ v h flv = -16 it m+T^, (3) 

where R^avp is the Riemann tensor of the MBH back- 
ground geometry. However, according to this equation, 
the metric perturbations diverge at the particle location. 
Then, the gravitational backreaction on the particle mo- 
tion, the self-force, is provided by the regularized part 
of the perturbations, say according to a Hadamard 
prescription given in |24| . Then, the equation of motion 
for an EMRI, the MiSaTaQuWa equation, is 

^+r>^ = m7/ 1 F«, (4) 
where the self- force, F" F , is given by 

F s a F = ~\m* (g aX + u a u x ) u»u v (2VXa - V A /i*„) . 

(5) 

Therefore, the dynamics of EMRIs is determined by the 
coupled system of Eqs. Q,([5|, and ([2| with a practical 
regularization scheme (like the mode sum scheme [75]). A 
remarkable point is that Eqs. Q and |5]) can be rewritten 
as geodesic equations of motion for a point particle in 
a perturbed geometry that only take into account the 
regularized part of the metric perturbations, i.e. geodesic 
in the geometry g af3 + h* p [75]. 



But how do we evolve the trajectory of the SCO ac- 
counting for the self-force in a self-consistent way? As it 
was proposed in [55J (see also [77] and [78] for a recent 
use of this technique) , one can use a relativistic extension 
of the well-known method of osculating orbits. The idea 
is to consider always the geodesic tangent to the trajec- 
tory, so that the motion transitions from one geodesic to 
the next. Such smooth transition is facilitated in EMRI 
by the fact that EMRI trajectories are very close to a (lo- 
cal) geodesic for a long time. This is because of the clean 
separation in EMRI time scales: the radiation-reaction 
time scale is much larger than the orbital (geodesic) one, 
except for the tiny fraction corresponding to the merger- 
plunge phase. 

The way to carry out this transition is to properly ac- 
count for the time-evolution of the set of orbital elements 
that completely characterize a geodesic orbits. Follow- 
ing [65] . it is important to distinguish between two sets 
of orbital elements: principal orbital elements, in our 
case either X A = (E, L Z ,Q/C) or O a = (e,p, t/6» inc ); and 
positional orbital elements that determine the initial po- 
sition in the geodesic as well as the geodesic initial spa- 
tial orientation. The radiation-reaction changes in the 
principal elements are due to the dissipative part of the 
self-force, while radiation-reaction changes in the posi- 
tional elements are due to the conservative part of the 
self-force. In this work we only consider dissipative ef- 
fects and hence we are only concerned with changes in 
the principal orbital elements. 

The implementation of the method of osculating or- 
bits consists in the translation of the fact that at any 
time there will be a geodesic trajectory, z£, with orbital 
elements whose position and velocity at that time will 
coincide with those of the accelerated trajectory. This 
can be written in the following way: 

z%t) = z«{t;V a {t);1 A {t)) , (6) 

^(t) = ^(t;V a (t);Z a (t)) , (7) 

where V A denote the positional orbital elements. Al- 
though we have used here the constants of motion as prin- 
cipal orbital elements, we could also have used O a . Com- 
bining these osculation conditions with the equations of 
motion Q , we can arrive at the following equations [55] : 

dz% dV A <9zg dl A 

cYP A ^ + &I A ^ ~ ' () 

d (dz°\ dV A d (dz%\ dl A _ 
dV A \ dr J dr + dl A \dr J dr SF ' 1 ' 

where a" F is the SCO self-acceleration, which is related to 
the self- force of Eq. ^ by a" F = m~ 1 F° F . This is due to 
the fact that X A has been defined as the SCO constants 
of motion per unit mass. From the inversion of these 
equations we can obtain the evolution of the different or- 
bital elements. In this paper, we will focus exclusively 
on the dissipative effects of the self-force which only af- 
fect to the principal orbital elements (either X A or O a ), 
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i.e. dV A /dr = and we ignore the first term in Eqs. Q 
and ^}J. We will refer to the principal orbital elements 
simply as orbital elements. 



III. THE NEW KLUDGE SCHEME 

In this section we present all the details of the new 
kludge scheme. As explained in the previous sections, 
we separate the problem into two parts: that of con- 
structing the SCO's trajectory and that of building the 
waveform from these trajectories. We begin then with 
a description of the background MBH geometry and its 
associated geodesies. We then show how to enhance the 
geodesic equation system to allow for the variation of the 
constants of the motion. The latter require knowledge of 
the self-acceleration, which we calculate in a multipo- 
lar, post-Minkowskian expansion. Finally, we present an 
explicit transformation between Boyer-Lindquist coordi- 
nates (used to evolve the modified geodesic system) and 
harmonic coordinates, needed to generate waveforms in 
a multipolar decomposition. 



where k a and l a are the two null principal directions of 
the Kerr geometry 



k a = 



A 



r = 



A 



1,0 



n 
A 



(13) 

When the effect of the self-force is neglected (or equiv- 
alently, in the limit of zero mass for the SCO), the 
SCO follows geodesies orbits of the MBH geometry. The 
Boyer-Lindquist coordinates of a timelike geodesic can 
be parameterized in terms of proper time as z^{t) = 
(t(r),r(T),6(T), </>(t)). For a given geodesic, we can con- 
struct three geodesic constants of motion, corresponding 
to each of the three symmetries of the Kerr spacetime. 
Stationarity leads to a conserved energy, E, or equiva- 
lently an energy per unit mass 



E = E/m* = -Ci'V 



(14) 



Axial symmetry leads to a conserved component of the 
angular momentum vector (the one along the spin axis, 
which we choose to be the z axis) 



(15) 



A. MBH Geometry and Geodesic Motion 

The geometry of the MBH is modeled by the Kerr met- 
ric |79| . a vacuum stationary and axisymmetric space- 
time that describes the final state of gravitational col- 
lapse, according to the BH no-hair conjecture |80j and 
uniqueness theorems (see, e.g. |81j). In Boyer-Lindquist 
coordinates [82], (zbl) = (t, r , 9, 0), the line element cor- 
responding to the Kerr metric, g^«, is given by 



ds 2 = -dt 2 + ^-dr 2 + p 2 d6 2 + (r 2 + a 2 ) sin 2 i 



asm 



(10) 



where p 2 



a 2 cos 2 I 



and A 



2M.r + a 2 



r f + a , with / = 1 — 2M,/r. For convenience, we also 



define the quantity £ 2 = (r 2 
function A has two roots: 



a 2 ) 2 - a 2 A sin 2 . The 



The Killing tensor symmetry (12 1 leads to a conserved 



Carter constant, which can be defined as follows 

Q = Q/ml = C ap u a u? , (16) 



but also we will use the alternative definition 



C = Q - (L z - aEY 



(17) 



The existence of these three symmetries makes the 
geodesic equations for z^(t) integrable and completely 
separable |83j . The separation can be carried out using 
the definitions of these constants of motion and also the 
normalization condition of the four-velocity, g^u^u 1 ' = 
— 1. The separated equations of motion for the compo- 
nents of z^(t) obey the following set of ordinary differ- 
ential equations (see [H] for a detailed analysis of the 
physical properties of the solutions of these equations) : 



p 2 — = \ (Y?E - 2M.arL 7 
dr A v z 



(18) 



M. ± v/m|" 



(11) 



The root r + (> r_) coincides with the location of the 
event horizon. 

The Kerr geometry is stationary, as described by the 
timelike Killing vector field £2 x = 8f, and axisymmetric, 
as described by a spacelike Killing vector field £?U = 8^ . 
It also well-known that the Kerr geometry has an ad- 
ditional symmetry described by a 2-rank Killing tensor, 
£ a p, given by 



to 







Ak (Jf3) 



(12) 



dr 



= [(' 



dOY 



E-aL z ] 2 - (Q + r 2 ) A, (19) 



1 

A 



C - cot 2 OLl 



2M. a rE 



2 2 

a cos i 



l-E 2 



sin 2 6 



A 



a 2 sin 2 i 



(20) 



(21) 



We are interested here in timelike bound and stable 
geodesies, what we call orbits. These orbits, apart from 
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being characterized by the three constants of motion I A 1 
can also be characterized by the orbital elements O a , 
which can be defined in terms of the turning points of 
the radial and polar motion (see Appendix [E] for more 
details). The turning points for the radial motion are 
just the minimum and maximum of r, also known as 
pericenter and apocenter, r por . and r apo , which can be 
used to introduce the concepts of semilatus rectum and 
eccentricity 



pM. 

1 + e 



pM. 

1-e 



or equivalcntly 
2r 



MJr . + r ) ' 

• \ peri 1 apo / 



r 

P' 

■ 7\ 



(22) 



(23) 



The turning point for the polar motion is just the min- 
imum of 9, min 6 [0,7r/2], which determines the inter- 
val in which 9 oscillates, i.e. (6* min ,7r — 6* min ) and it can 
be used to introduce the concept of inclination angle, 
° in a e [-7r/2,7r/2], as 



0i„c = sign(i 2 ) 



(24) 



Another common definition of orbital inclination angle 
uses the constants of motion X A 



COS L — 



(25) 



Alternatively, the orbits can also be characterized in 
terms of three fundamental frequencies (see e.g. [851 - 
[87] ) with respect to the Boyer-Lindquist coordinate time 
(they can be also constructed using proper time or any 
other time): f2 r , associated with the radial motion (from 
periapsis to apoapsis and back); f2 e , associated with po- 
lar motion; and £2^, associated with azimuthal motion. 
These frequencies are important because precessional or- 
bital effects are due to mismatches between them and 
because they can be used to decompose, among other 
things, the GW in a Fourier expansion [86] . 

Expressions for these frequencies in terms of quadra- 
tures have been obtained for Kerr in [85] . and recently 
also in [86] [87] . Appendix [F] provides explicit formulas 
for the fundamental frequencies and periods used in this 
paper. 

A final consideration regarding geodesic motion that 
is going to be important in this paper is the 3 + 1 split- 
ting of the four-velocity into spatial and time compo- 
nents. This splitting is associated with the time vari- 
able used to evolve the geodesic equations, in this paper 
Boyer-Lindquist time t. The four-velocity of the SCO, 
u a = dz^/dr, choosing the Boyer-Lindquist time param- 
eterization z^(t) — (t,z' l (t)), can then be decomposed as 



dt dz l 
dr dr 



= r(Ly) : 



(26) 



where T and v % are given by 

r = ««=* ^ = p-V = ^ 

dr dt 

Using the normalization condition gjj^u^u 1 ' - 
factor r, the GR generalization of the special-relativistic 
Lorentz factor, can be written in terms of the metric and 
the components of the velocity v % as follows: 



(27) 



-1, the 



F = (-g«-2gfy- g K 



-1/2 



(28) 



B. New kludge Osculating Trajectories 

In this work we consider a version of our scheme in 
which we only include the dissipative effects of the self- 
force, i.e. those that only affect to principal orbital ele- 
ments. In Sec. [VI] we discuss how to introduce conser- 
vative pieces of the self-force in the new kludge scheme. 
The time scale of change of X a (t)/O a (t), the radiation- 
reaction time scale T nni is much bigger than the orbital 
time scales T_ , .. , , such that the ratio of these time scales 

Orbital ' 

satisfies: T 0rbital /T RR - q. 

Following the method of osculating orbits described in 
Sec. [n] we can describe the orbital evolution as given 
locally in time by the geodesic Eqs. (18)-(21), where 



the constants of motion are promoted to time-dependent 
quantities: 



Z A (t), 



O a (t) 



(29) 



Although here we parameterize the time dependence 
of these quantities in terms of proper time, in prac- 
tice we will use coordinate time, that corresponds to 
time associated with distant observers. For our discus- 
sion, we denote the solution of the evolution Eqs. ( 18 1- 



(211 with the modifications of Eq. (29 1 by z^(r) 



[t(r), r(r), 9(t), 4>(t)], but it is clear that it would no 
longer be a solution of the original geodesic equations. 
We do not decompose this solution, which takes into ac- 
count self-force effects, into a background geodesic plus 
a deviation |26j . as the deviations grow secularly in time 
and after a certain number of cycles it cannot be consid- 
ered a small deviation of the background geodesic orbit. 
Instead, in the spirit of the osculating orbits method, we 
treat z^ as a new, self-consistent trajectory that is contin- 
uously corrected away from geodesic motion during evo- 
lution. In this sense, we are constructing an orbit that is 
made out of geodesic patches corresponding to different 
constants of motion. The transition from one patch to an- 
other is given by the (multipolar, post-Minkowskian) self- 
force, lacking a more accurate prescription. The length 
(duration) of the geodesic patches, or equivalcntly, the 
frequency at which the constants of motion are updated, 
is a free parameter that we can choose in the new kludge 
numerical implementation. 

The evolution equations for X a (t)/O a (t) can be ob- 
tained from the equations of osculating orbits, Eqs. ([8l 
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-C [t) a a 

Sa u sf ' 




(30) 


c(<t>) a a 

Sa "sf ' 




(31) 


2 £af3 U<,£t SF l 




(32) 


^+2(o£?-L.) 


\ dr dr J 


(33) 



and (J9j) . In our case, the inversion of these equations to 
find dX A (r)/dT can be done easily by using the symme- 
tries of the Kerr geometry. By applying g^Cf t) , SapC^y 
and ^pU 13 to Eq. ^ in combination with Eqs. (|l~4|)-(fl7|) 
we obtain the evolution equations for E, L z , and C/Q 
respectively. 

dE 

dL* 
dr 
dQ 
dr 
dC 

The evolution of these quantities with respect to coordi- 
nate time t introduces a factor in the above equa- 
tions, due to the relation: d/dt = T~ 1 d/dr. The most im- 
portant quantity here is the self-acceleration a" F , which 
as we shall see scales ~ T 2 [see Eq. ( [57] )] and is given by 
the metric perturbations h^p. This is probably the main 
ingredient of self- force descriptions of EMRIs, which is 
described in more detail in the next section. 

For certain special orbits, the rate of change of the or- 
bital elements satisfies special relations due to the sym- 
metries. One of these are equatorial orbits, i.e. orbits 
with 0(t) = n/2. We can see that from Eq. (20) this 
implies C = 0. Therefore, equatorial orbits are charac- 
terized by 6(t) — ir/2 and C = 0. Equation (33) allows 
us to write dC/dr — 2C a af F , where 



{€ + u a u' 3 ) [^ x u X + (aE - Lz) + 



(34) 

The vector C a vanishes for equatorial geodesies, which 
implies that the Carter constant C(t) is always zero and 
Q(t) can be obtained directly from a combination of 
dE/dt and dLz/dt. The self-force then maps equatorial 
geodesies to equatorial geodesies, i.e. the equatorial char- 
acter of geodesies is preserved upon self-force evolution. 

Another type of orbits with extra symmetries are equa- 
torial and circular orbits. These orbits have a helical 
symmetry described by an approximate Killing vector 
(which is exact on the geodesic intervals). This symme- 
try can be derived from the fact that for these special 
orbits d(j)/dt = fL = const, [see Eq. (21)]. The angular 
velocity fL can be written as (see [84]): 



a, = ±- 



M. 



1 ± — v 

1 ^ M. u ( 



(35) 



where the upper/lower sign corresponds to pro- 
grade/retrograde circular orbits (i.e. that coro- 
tate/counterrotate with the MBH spin and have 

L z > ®/ L z < 0): v o = V M »/ r o, and r o is 
the Boyer-Lindquist radial coordinate of the cir- 
cular orbit. Then, the helical Killing vector is 



C, = C ( Q t) + fyMCfo = d? + %{r)d%. The associated 
constant of motion is then A = g^CHci -1 ^ = — E + ^lu > L z . 
The evolution of A is dA/dr = Ca ela sF an d it can be 
shown that dA/dr = for locally circular-equatorial 
orbits. Hence, the evolution of the angular momentum 
in the spin direction is related to the evolution of the 
energy by dL z /dt = il^ 1 dE/dt. 

Finally, let us consider circular but nonequatorial or- 
bits, i.e. orbits with r = r a = const, but C ^ 0. There is 
no helical symmetry for these orbits due to the fact that 
the MBH spin is not aligned with the orbital angular mo- 
mentum. Nevertheless, as shown in [55JIHS] the radiation- 
reaction evolution of these orbits preserves their circular 
character (see also [42]). These orbits are then character- 
ized by the vanishing of the right-hand side of Eq ( |19| ), 
R(r ) ee [{rl + a 2 )E~ aL z ] 2 - (Q + r 2 ) A(rJ = (the 
radial coordinate does not change), and its first radial 
derivative, R'(r ) = (dR(r)/dr) r = (the orbit is al- 
ways at a turning point of the radial motion). In ad- 
dition, the condition R"(r a ) < has to be satisfied for 
the circular orbit to be stable. Following [J3] , the preser- 
vation of the circularity of the orbit along the inspiral 
translates into the following two conditions: R — and 
Rf = 0. From these two conditions we can express the 
time evolution of the radius of the orbit and the evolu- 
tion of the Carter constant in terms of the evolution of 
the energy and angular momentum of the orbit: 



-21 °22 



E 
L, 



(36) 



where the coefficients c AB (A,B = 1,2) and d are func- 
tions of (M ml a; E, L z , C, r a ) and are given in Appendix [G] 
Therefore, the evolution of C and r a can be obtained by 
evaluating dE/dt and dL z /dt from the multipolar, post- 
Minkowskian self-force. 

To finish this section, and as a consistency check, we 
take the Newtonian limit of the rate of change of the 
orbital elements in Eqs. (30)-(32). In this limit, we find 
that to leading order 



771 Z I 



r 2 sin 2 a 



Q = -2r 2 sin 2 



C 



RR 3 

a a RR , 
AM.ar cos 2 6 sin 2 9 a* 



(37) 
(38) 
(39) 
(40) 



which agrees with the Newtonian expressions of [90] after 
transforming to Cartesian coordinates. Notice, however, 
that Eqs. ( 30 )-( 32 ) contain many more terms than in [SU] , 
due to the fact that in the latter the T factor is expanded 
in the low-speed and weak-field limit. 



C. Multipolar Post-Minkowskian Self- Acceleration 

Here we discuss a method to obtain an approxima- 
tion for the self-force The most rigorous approach 
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would be to compute this force within the framework 
of BH perturbation theory, but as already argued, such 
a task is still under development and in computational 
terms would be very costly. A different approach is to 
extract the self-force from the PN equations of motion. 
Such a path was taken by Pound and Poisson [65] for 
a Schwarzschild background and Gair, et al. for a Kerr 
background [78j . This approach, however, is built under 
a PN approximation, which is not ideal for EMRI mod- 
eling as the SCO can reach relativistic velocities, with 
significant relativistic T factor, as it orbits close to the 
MBH. 

Instead of either of these approaches, we here approx- 
imate the self-force via a multipolar, post-Minkowskian 
expansion (see e.g. [6"4"l |9"TI |9"2"] ) : 



PM 

=00 



PM 



-\ + 2V -2V 2 + 2V T 
-AV* - W* R + 0(G 



p p 

9/2 



0(G 9/2 



). 

g™ = 6 lJ [l + 2V + 2V BIi +0(G 4 



(41) 
(42) 
(43) 



where V and V 1 are time-symmetric potentials, in the 
sense that they are made out of the half-sum of re- 
tarded and advanced integrals of the stress-energy ten- 
sor over the source, i.e. the half-sum of the retarded 
and advanced Green functions associated with the per- 
turbative equations. As a consequence, these potentials 
are invariant under time inversion. The quantities V nn 
and V£ R are time- asymmetric radiation-reaction poten- 
tials, constructed from the half-difference of retarded and 
advanced waves, and hence, odd under time inversion. 
These potentials are given by [BH [93] 



189 



~xlx%Ml]\t n ) + 0{xlM\f ), (44) 



VL(t*,*„) = ~ x n jk>M % ) ( t «) 



21 

4 

45 



..x^S^it^+OixTM^), (45) 



where we recall that (x^) — (i H ,a4) are spacetime har- 



monic coordinates, 
symbol, and 



£j -j. is the antisymmetric Levi-Civita 



~<ijk> _ ijk 



„2£(y fc) 



(46) 



is the symmetric trace-free (STF) projection of the multi- 



index quantity x 



ijk 



X X J X 



ix k . The first term in V„ 



[Eq. ( 44 )] corresponds to the well-known Burke-Thorne 
radiation-reaction potential 67J: 



certain derivatives of the asymmetric sum (half-difference 
of retarded and advanced waves) of multipole moments 
(see, e.g. Eqs. (2.8) in [53] )■ Equation (45), however, 
is obtained by expanding these integrals in a slow- 
velocity approximation, after which the arguments of the 
radiation-reaction potentials depend on time only. This 
is consistent with the fact that the radiation-reaction 
force is to be evaluated in the source zone of the SCO, 
and not in the wave zone. 

Let us provide explicit expressions for these multipole 
moments. To lowest order, the mass moments are given 
by 



M s , 



hj ~ r l m Z <ij> ' lvJ -ijk 

and the current moment is 



rj 8m z 



<ijk> 



Sij — rj5m e kl<i Zj > k z l . 



(48) 



(49) 



where angle-brackets are STF projections. To higher or- 
der, these moments become more complicated as there 
are now nonlinear contributions from the nonreactive 
potentials (i.e. nonlinear contributions from the back- 
ground) as well as tail and memory contributions. These 
expressions can be found for example in Eq. (3.1)-(3.3) 
of [94] and Eq. (5.3)-(5.5) of [95]. In BH perturbation 
theory language, these higher-order terms would con- 
tribute conservative and dissipative corrections to the 
dissipative equations of motion. We neglect these con- 
tributions in the current version of the new kludge ap- 
proach, but they can be easily incorporated in future 
improvements of the scheme. 



The metric given in Eqs. ( 41 )-( 42 ) is expanded in the 



far-field limit, a resummation of which is necessary in or- 
der to use it for self-force calculations. All terms that are 
independent of the radiation-reaction potentials can be 
identified with MBH background geometry terms, cor- 
rected perhaps by the presence of the SCO. Neglecting 
the latter, we can exactly resum the metric in Eqs. (41)- 
( 43 1 so that it can written as follows 



PM 
itt 

J.PM 
J.PM 



K,H 
Itt 
,K,H 

oti 

,K,H 

lij 



Kt 



G(G 9 / 2 ), 
G(G 9 / 2 ), 
0(G 4 ), 



(50) 
(51) 
(52) 



where g^ H is the Kerr metric in harmonic coordinates 
(we shall discuss the issue of coordinates in more detail 
in Sec. HID) and where we have introduced the metric 



perturbation, hf^, which is given by 



. RR 

n tt 



2K 



. RR 

n ti 



-4V 1 

1 * T3 ■ 



2^Kr 



(53) 



Burke— Thornc 



(t s , X H ) 



1 



The quantities M^ n) , M$, and S£> are the nth- 
time-derivative of the STF mass quadrupole, mass oc- 
topole and current quadrupole moments. Formally, the 
radiation-reaction potentials depend on the integral of 



(n) 



(47) 



These are the metric perturbations that are going to de- 
scribe the radiation-reaction effects in our new kludge 
scheme, and hence they are our approximation to the 
regularized metric perturbations h*g in Eqs. ^ and (|5j). 

Regarding the MBH background geometry, it is useful 
for the purposes of this work to introduce some scalar, 
vector, and tensor potentials in harmonic coordinates, 
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which contain information equivalent to the potentials 



V and V 1 in Eqs. (41)-(43) 



K = K 00 =g«f + l, Q^Q°° = g°° H + l, (54) 



t — S>0i i 



<T = g£ H . (55) 



K = o- K,H 
ij = Sij 



Q 



n = 



>K,H 



(56) 



or more compactly i£ = g^ H — ?7 and Q My = g 



/I /.' 



The expressions from these local potentials can be 



derived directly from the expressions of the components 
of the Kerr metric in harmonic coordinates given in Ap- 
pendix |B2j 

We are now in a position to compute the self- 
acceleration in terms of the radiation-reaction potentials 
of Eqs. ([44]) and |45} through Eq. ([53]) and Eq. (5). In 
what follows we develop the expression for the acceler- 
ation in order to describe its structure and the role of 
the different terms that appear on it. The first step is to 
use the decomposition of the SCO four-velocity given in 
Eqs. (26) and (27) in such a way as to rewrite Eq. ^ in 



the following form: 



= -FT' 



aB 



(57) 



where we have factored out the dependence on the rela- 
tivistic r factor [see Eq. ( A3 ) for the expression of V in 
terms of the potentials K„ v and the spatial velocity v % \ 
P al3 is the projector orthogonal to the SCO four-velocity, 
pa/3 _ g«/3 _|_ yayfi ^ m terms of the potentials Q^ v it is 
given in Eq. ( A2 )) . The two pieces of the self- acceleration 



in Eq. (57) contain different terms. The first piece, A a 



(i) 



only contains gradients of the radiation-reaction poten- 
tials and the spatial velocity v % , whereas the second piece, 

(2) 

A a L contains explicit couplings between the MBH po- 
tentials (actually their derivatives) and the radiation- 
reaction potentials. The form of these two terms is the 
following: 



4(1) = cm v f v " 



where 



and 



A (2) _ _J,BH p/3 a v 



RR\ 



(58) 



(59) 



(60) 



where the connection is to be computed with the Kerr 
background metric in harmonic coordinates. Putting all 
these different ingredients together and separating space 
and time components, we arrive at the following structure 



for the piece A. 



(i). 



A (l) = A KU 



A 1 ) _ A RR. 



(61) 



where the expressions for the quantities A RR and Af R , 
which depend only on the three- velocity v l and spacetime 
derivatives of the radiation-reaction potentials V RK and 
V* R , are given in Eqs. (A4) and (A5) of Appendix [A] 



The second piece of the self-acceleration, A 
following structure: 



(2) 



V-RR 
i 



Af> = -2 [B RR + C RR + V RR ] V- 



2 5 * 



RR 



c 3 

RR 



V J V 

RRl v RR ! 



has the 



(62) 



(63) 



where the quantities B RR , B RR , C RR , C RR , £> RR , and £> RR , 
which depend on the SCO three- velocity v l and the MBH 
potentials (K, K i: K^) and (Q,Q l , Q %] ) and their spatial 
derivatives, are all given explicitly in Eqs. (A6)-(All) of 
Appendix [X] 

One can check that in the limit gjjjf — > r/^ (or equiv- 
alently, (K, K i ,K lJ ) -> (2^(1-^,-4^,2^) and 
\v\ <C 1 (in units c — 1), the spatial components of 
the acceleration in Eq. ( [57] ) reduce exactly to Eq. (3.11) 
in [911 [92] order by order. This can be checked by realiz- 



ing that the self-acceleration here, Eq. (57), is related to 
the self-acceleration in [211 [§2] > let us call it a\ 



via 



dt 2 



~dt 



= r- 2 (4 



,)• 



(64) 



D. From Boyer-Lindquist to Harmonic Coordinates 



The resummed MBH background metric in Eq. (52), 



the background Kerr metric, must be written in harmonic 
coordinates in order for all terms to be in the same co- 
ordinate system and to be consistent with other ingredi- 
ents of the new kludge scheme, like the waveform gener- 
ation procedure described in Sec. |IIIE[ However, there 
are parts of this scheme that are easier to implement in 
Boyer-Lindquist coordinates, such as the integration of 
the geodesic equations. Therefore, it is crucial to de- 
termine out how to transform from Boyer-Lindquist to 
harmonic coordinates. 

Harmonic coordinates refers to any member of the fam- 
ily of coordinate systems, say {a;^}, that satisfy the equa- 
tion 



0. 



(65) 



where □ is the D'Alembertian operator: □ = g Q ^V Q V /3 . 
When we write down the D'Alembertian operator in this 
coordinate system, this condition is equivalent to the 
condition: E^^aB = ^' which m turn is equivalent to 

requiring <9 (3 Sh' 3 = 0, where = ^/-g H Eu ■ Notice 
that the harmonic coordinate condition and the harmonic 
gauge condition are different things. The first expresses 
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a property of a given coordinate system in a given space- 
time, whereas the second refers to the correspondence 
between the background and perturbed spacetimes in 
perturbation theory. In this sense, the gauge condition 
d^h^ = 0, where ft/ 1 " = g a/3 — rf 1 ^ , enforces harmonic 
coordinates at first order in perturbations around flat 
spacetime. Although the Lorenz gauge condition Eq. ^ 
is similar to the harmonic gauge condition, these condi- 
tions are in general different, coinciding only for pertur- 
bations around Minkowski spacetime. 

The most commonly employed coordinate systems in 
GR to describe BH geometries are Schwarzschild or 
Boyer-Lindquist coordinates, neither of which is actu- 
ally harmonic. Although the time and azimuthal Boyer- 
Lindquist coordinates, t and <p, satisfy Eq. ( 65 ) , the radial 



and polar coordinates, r and 9, do not. One can further 
check that Eddington-Finkclstcin and Kerr-Schild coor- 
dinates are also nonharmonic. Of course, if one consid- 



ers a vacuum, Minkowski spacetime, then Eq. (65) be- 
comes trivial, d a d a x fl — 0, which is generically satisfied 
by many coordinate systems. 

Harmonic coordinates for the Kerr metric have been 
studied extensively in the literature. A harmonic slicing 
of the Kerr metric was found in [96] , where the time- 
function T is related to Boyer-Lindquist coordinates via 



M. 



V2 



+ y(r 2 - a 2 ) + 4a 2 z 2 



1/2 



9 = arccos 



r - M. 



We have here introduced the short-hand notation 



r B =(xl + y*+4) 1 ' 2 .. 



and the angle function ^(r), which is given by 



W 2 ll-^ntr)/' 



with 



fi(r) = tan 



2^W 



In 



(72) 



(73) 



(74) 



(75) 



(76) 



T = t 



hi 



(66) 



where r ± are given in Eq. (111. Although the slicing is 



harmonic, the full set of four-dimensional coordinates is 
not. Ref. [HZ] did find a full set of harmonic coordinates 
for Kerr but no explicit expressions for the metric com- 
ponents were given. In [98}, based on work by Ding (see 
Ref. [95] for references on this), a different transforma- 
tion to harmonic coordinates was found, where the time 
slicing was not modified. This is very convenient for the 
new kludge scheme, which is why we adopted this choice 
in this paper. 

Following [98] . we can map the Kerr metric from 
Boyer-Lindquist coordinates (t, r, 9, <p) to harmonic ones 
(t u , x a , y u , z H ) via the coordinate transformation 



(67) 



The coordinate transformation between Boyer-Lindquist 
and harmonic coordinates requires expressions for 
(cos $, sin $) and the Jacobian and Hessian of the trans- 
formation, which we provide in Appendix |B| 

We have checked that Eq. (65) is satisfied identically 
for (x") = (t H , x H , y H , 2 H ) in the Kerr background. More- 
over, one can see that this transformation reduces to 
the standard one from Schwarzschild to harmonic coordi- 
nates [9 9 1 in the a — > limit. In this limit, the transfor- 



mation ( 68 )-( 70 ) reduces to the Euclidean transformation 
from spherical to Cartesian coordinates, with r H playing 
the role of the spherical coordinate, related to the radial 
coordinate r by r H = r — M % . 

This explicit coordinate transformation allows us to 
compute certain ingredients of the new kludge scheme in 
Boyer-Lindquist coordinates and then transform the re- 
sult to harmonic coordinates when needed. For instance, 
the Kerr metric in harmonic coordinate is then simply 



x H = J (r - M.) + a 2 sin cos [0 - $(r)] 



(68) 



r K,H 



r K,BL " A BI- ^-"BL 



(77) 



y n = y(r- M.) + a 2 sin 9 sin^ - $(r)] , (69) 



z H = (r — M,) cosf 
while the inverse transformation is 



t = t 



(j> = $(r) + arctan 



(70) 



(71) 



where the Jacobian dx^/dx^ is given in Appendix B 1 



and the transformation is shown explicitly in Ap- 
pendix |B2l 

In particular, in the current implementation of the new 
kludge scheme, we choose to perform the quasigeodesic 
evolution using variables associated with Boyer-Lindquist 
coordinates (see Sec. IV). Then, from the information 



of the trajectory in Boyer-Lindquist coordinates we can 
compute the velocity and accelerations in harmonic co- 
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ordinates using the following relations 



dxL 



dr 



dr 



dr 89 
d 2 r„ 



dr 

+ ?T^ (78) 



BL 1 r\ ^ r\j BL BL 



dr 
d 2 r 



W 6 

d 2 r 



BL BL 

dr„ 



d 2 9 



a 2 r H 

d 2 (f> 



d 2 r 

( d 2 T ■ d 2 v ■ d 2 v 

2 I rO+Z^r r<t> 



\ dr 80 



dr dtp 



d6d<j) 



where the Hessian d 2 r n /dx Blj dx BL is given in Ap- 
pendix |B 3| Notice that we only need the spatial com- 
ponents of the self-acceleration in Eqs. ( 30 )-( 32 ). With 



this acceleration, we can then compute the change in the 
constants of the motion as we osculate from one geodesic 
to the next. 

On the other hand, thinking about extending the new 
kludge scheme to systems in which the central object, the 
MBH in our case, is not described by the Kerr metric but 
by a different metric tensor (either because the central 
object corresponds to an exotic distribution of mass or 
because it is governed by an alternative theory of grav- 
ity, or perhaps both), it may happen that transforming to 
harmonic coordinates may be a very difficult task, even 
when given an explicit coordinate transformation as the 
one above is known (although in this case one may re- 
sort to numerical computations). Fortunately, one does 
not need a full coordinate transformation to read off the 
multipole moments of the background. Thorne [47] has 
shown that it suffices to work with so-called asymptoti- 
cally Cartesian and mass-centered coordinates of order N 
(ACMC-iV). These coordinates are defined such that the 
background metric is time-independent and has a certain 
spherical harmonic structure (see, e.g. Eq.(ll.la) in [47] ) . 
Thorne has further shown that harmonic coordinates are 
ACMC-oo, and found an explicit map between Boyer- 
Lindquist and ACMC-2 coordinates [17] . In Appendix [C| 
we construct an ACMC-4 coordinate system and com- 
pare it to the exact harmonic coordinates described in 
this section. The difference between a standard ACMC- 
4 and its harmonic generalization appears only near the 
horizon, and in particular, it should affect results when 
PN corrections to the multipole moments are taken into 
account. We work here directly in harmonic coordinates 



as given by Eqs. (67 1- (70 1 and explained in detail in Ap- 
pendix [5] 



E. Waveform Generation 

Once the orbital evolution, including self-force back- 
reaction, has been computed, one can construct the re- 
sulting GWs. Expressions for these as a function of 



the trajectories can be obtained by solving Eq. (|3| ei- 
ther numerically or analytically, via some approximation 
scheme. Here we use the solution for the waveforms that 
comes from the combination of post-Minkowskian and 
PN expansions and which gives rise to standard multipo- 
lar expressions (see e.g. |100j for a review). In this way we 
can write the plus- and cross-polarized solution to Eq. ([3]) 
as an infinite series expansion in terms of derivatives of 
mass and current multipole moments |47L 1 101 j : 



i ky 



(80) 



■ (79) 

where t r denotes retarded time and U lm and V 



Cm 



arc- 



radiative mass and current multipole moments: 



U 



16tt 



(2£ + 1)!! V 2£(£-l) 



(£ + !)(£ + 2) tm , 

u l y l > 



(81) 



V 



tm 



-32n£ 



(2£ + l)!!V 2£{£ + !){£-!) 



v L yr\ (82) 



where (Ul,Vl) are symmetric, trace- free, mass- type and 
current-type multipole moment tensors (see, e.g. Eq.(15) 
in |101j ). Even for the case of comparable-mass BH merg- 
ers, such an expansion taken only up to quadrupole order 
(/ = 2) has been shown to be sufficient to recover the 
waveform to excellent accuracy [102j . 

This prescription is then complete, once expressions for 
the radiative moments are given in terms of derivatives 
of the orbital trajectories. Such identification, however, 
is difficult, as these moments are defined in the far-zone 
(many gravitational wavelengths away from the center 
of mass of the binary), and have no knowledge of the 
source multipole moments, defined in the near-zone (less 
than a gravitational wavelength from the center of mass) . 
Asymptotic matching can be used to relate the radiative 
to the source multipole moments |100j . yielding explicit 
expressions that depend only on derivatives of the orbital 
trajectories in the near-zone. Here, in this initial version 
of the new kludge scheme, we only consider just the lead- 
ing order contributions to the multipoles moments, i.e. we 
identify the source and radiative moments: Lt L — M L and 
V L = S L , ignoring in this way subleading corrections that 
correspond to tail and memory effects. 

In a transverse-traceless gauge, one can rewrite the 
harmonically decomposed metric perturbation in the fol- 
lowing simpler form: 



oo 

£ 

1=2 



--M (t ' 
£\r L ~ 



i(t r )N 



L~2 



1 



t kl(iSj)kL- 



-l(*r 



)n t N L _ 2 



(£ + 1)! r" M " ' L ~ 1 ' ' ' ' L ~- ' ^ 

which, when we consider only multipoles up to the mass 
hexadecapole and the current octopole, reduces to 
2 



2 ■■ 

-M STF 

l 

Gr 



3r 



Mijkn" +Ae kl {i S ])k ni 



[Mijkin k n l 



6e kl {t Sj)kmn l n r ' 



(84) 
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The expressions for these multipole moments have been 



given in Eqs. (48) and (49), except for the mass hexade- 



capole and current octopole multipoles, which are given 
respectively by 



M, 



ijkl — V m z <ijkl> ' 



5' 



ijk — V 771 e lm<i Z jk> 



I •m 



(85) 



(86) 



Again, these moments contain higher-order corrections 
that can also be easily included in subsequent improve- 
ments of the new kludge scheme. The plus and cross- 
polarized projections can then be constructed via 



,n i,tt 



(87) 



where 



is the plus- and cross-polarization tensors. 



Finally, the observables that one wishes to compute are 
the GW response functions, which are given by a projec- 
tion of the plus- and cross-polarized waveform with the 
beam-pattern functions of the detectors. For a detec- 
tor like LISA there are two such functions, F + j/h and 
Fx,i/ii (see, e.g. [TB I 1103) for expressions of these func- 
tions) and the response is thus 



V3 
2D T 



A=I,II 



F x,A h > 



where D L is the luminosity distance from the source to 
the observer and the prefactor of V3/2 is due to the tri- 
angular arrangement of the LISA detector. 



in Eqs. (C5)-(C8)] and reactive potentials, in turn given 
in Eq. (45). To evolve osculating orbits, one must there- 



fore map the Boyer-Lindquist acceleration and velocity to 
harmonic coordinates via Eqs. ( 78 )-( 79 ), so as to compute 



the rates of change of the orbital elements, which then 
in turn allows us to map between osculating geodesies. 
Once the SCO's worldline has been completely obtained, 
we can use it in harmonic coordinates in the waveform 
prescription. 

The kludge nature of the approach then becomes clear. 
We employ a combination of approximation schemes that 
include a multipolar, post-Minkowskian expansion (for 
the far-zone metric perturbation and for the local pre- 
scription of the self-force) a post-Newtonian expansion 
(for the multipole moments in terms of the trajectories) 
and a BH perturbation theory expansion (when treating 
the trajectories as self-adjusting geodesies). All of this is 
tied together via a nontrivial numerical implementation 
that is described next. 



IV. NUMERICAL IMPLEMENTATION OF THE 
NEW KLUDGE SCHEME 

In this section we provide details of how we have imple- 
mented each of the ingredients of the new kludge scheme 
described in the previous section. The numerical code 
that we have developed is written in C language |104j and 
uses different parts of the GNU scientific library |105j as 
we describe below. 



A. Integration of the Equations of Motion 



F. Summary of the New Kludge Approach 

The physical quantity one is interested in is the re- 
sponse function, which is given by Eq. ( 88 ) in terms 
of the beam-pattern functions (see, e.g. [TS]) and the 
plus- and cross-polarized waveform. The latter is given 
in Eq. (87) in terms of the transverse-traceless metric 



perturbation. The first approximation we make is to ex- 
pand hjf in post-Minkowskian multipole moments Mj», 



A'hjk and Sij through Eq. (84), where we neglect tail 



and memory corrections. The second approximation we 
make is to treat these moments in a Newtonian-like fash- 



ion through Eq. (48) and (49) in terms of the trajectory 



of the bodies, neglecting post-Newtonian corrections due 
to nonlinearities. These two approximations provide the 
response function as a function of the trajectory of the 
bodies in harmonic coordinates. 

The orbital trajectories are obtained by solving the 
geodesic equations enhanced by time-varying orbital el- 
ements in Boyer-Lindquist coordinates. The time varia- 
tion of the orbital elements is prescribed by the radiation- 
reaction acceleration in Eqs. (61), (62) and (63) in har- 



monic coordinates and in terms of quantities that depend 
on the harmonic Kerr background [with the map given 



We need to integrate numerically the set of ODEs con- 
sisting of Eqs. ( 18 )-(21 ), and at the same time we need to 



update the value of the constants of motion (E, L z , C/Q). 
The separation of the geodesic equations has produced 
equations for the radial and polar Boyer-Lindquist coor- 



dinates, Eqs. (19) and (20) respectively, that have turn- 



ing points (at pericenter and apocenter in the case of the 
radial coordinate, and at the location of the orbital incli- 
nation angle in the case of the polar coordinate) , and at 
these points we have either f = or 6 = 0. This means 
that these are not the best variables for the numerical in- 
tegration, as ODE solvers present convergence problems 
at turning points. 

To avoid this problem, we introduce new variables in 
the place of the radial and polar Boyer-Lindquist coordi- 
nates, r and 9. These new coordinates are angle variables 
defined by the following expressions: 



pM. 



1 + e cos ip 



cos 2 9 = cos 2 6> min cos 2 x • (89) 



We can write the equations for r and 9 in terms of their 
turning points and other extrema (points at which the 
time derivatives vanish but are not accessible to the mo- 
tion). In the case of the radial motion, we can write the 



16 



right-hand side of Eq. ( 19 ) as 

(1 - £ 2 )(r apo - r)(r - r pc J(r - r 3 )(r - r 4 ) , (90) 

where r 3 and r 4 satisfy r apo > r pcr . > r 3 > r 4 . In the same 
way, we can write the right-hand side of the equation for 
the polar motion [Eq. |20|)] in the following form 



where 



a 2 {l-E 2 ) 



1-z 



z = cos 9 . 



cos 2 (9 



(91) 



(92) 



and z + > z_. We describe in Appendix [E| the rela- 



tions between these extrema, i.e. (r a 
(z_, z + ), and how to find them. 



, Vr.' r 3; r 4) and 



For convenience, we parametrize the trajectory in 
terms of the Boyer-Lindquist coordinate time t, which 
is also a time harmonic coordinate, instead of the proper 



time r. This is done using Eq. (18) to rewrite the evolu 



tion equations with respect to t. In this way, the resulting 
equations for (^>(i), x(t), </>(t)) are (see also, e.g. [86]): 



dtp _ Af.%/1 - E 2 y/\p(l - e) -p 3 (l + ecosV')] b(l + e) -p 4 (l + ecos-0)] 
~dt = (l-e 2 )[$W+a 2 £z(x)] 



(93) 



d x ^/a 2 (l-£ 2 ) (z + -z_ cos 2 x) 
~dt ~ ■$>(il)) + a 2 Ez(x) 



(94) 



d<t> _ 1 J 2M.ar{i))E 

lit ~ ^(ifj)+a 2 Ez(x) { A(r(V>)) 



1 _ a 2 
l-z(x) A(r(V)) 



(95) 



where we have introduced the following definitions: p 3 
r 3 (l - e)/M„ p 4 = r 4 (l + e)/M„ 



*(V0 



„2\2 



A(rty>)) 



2M.ar(ip)L z 



(96) 



and r("0) and z(x) are given through Eq. (89). There- 



fore, the actual outcome of the numerical integration of 
the ODEs of Eqs. (93 )-(95 ) is a time series of the three an- 
gles (t(>(t), x(t)i 4>(b))i which grows monotonically in time. 
The numerical method we use to integrate these ODEs 
is the Bulirsch-Stoer extrapolation method ( |106j ) as de- 
scribed by [M] (see also [TPS] ). 



B. Estimation of Time Derivatives 

Probably the main challenge in the numerical imple- 
mentation of the new kludge scheme is the evaluation of 
the time derivatives of the different quantities involved. 
To understand the nature of this problem let us focus 
on the computation of our multipolar, post-Minkowskian 
self-force. If we look at the expression for the radiation- 
reaction potentials V RR and V^ R , Eqs. (44) and (45), 



we realize that we need to compute up to the seventh 
time derivative of the mass quadrupole and octopole mo- 
ments and up to the fifth time derivative of the current 
quadrupole moment. But since we also need to compute 



time derivatives of these potentials [see Eqs. ( 57 ) and ( 61 ) 



r 



and also Eqs. (|A4j) and (|A5j)] we also need the eighth time 
derivative of the mass quadrupole and octopole moments 
and the sixth time derivative of the current quadrupole 
moment. 

In principle, one could think about computing these 
derivatives analytically by using the equations of mo- 
tion, Eqs. ( 18 )-( 21 ) . The problem is that we need the 
time derivatives of the trajectory in harmonic coordi- 
nates, and to pass from the ODE angles (ip(t), x(i), <fi(t)) 
to harmonic coordinates we need to first use Eq. ( 89 ) to 



go from these angles to Boyer-Lindquist coordinates, and 
then Eqs. ( 67 1-( 70 ) to go from Boyer-Lindquist to har- 
monic coordinates. Therefore, we would need to obtain 
analytically higher time-derivatives of the ODE angles 
(up to eighth order) , which involves using the Christoffel 
symbols of the Kerr metric and several of their deriva- 
tives, and also to differentiate several times Eqs. (89) 



and (67)- (70). In practice, this makes the analytical com- 



putations unfeasible, even using modern computer alge- 
bra systems, and even if they were not, one would worry 
about the reliability of the numerical evaluation of the 
huge expressions that one would result. 

For these reasons, we resort to a numerical evaluation 
of these derivatives. The starting point is the fact that we 
can compute the trajectory, the velocity, and the accel- 
eration almost directly from the integration of the ODE 
Eqs. ( 18 )-( 21 ) and with a high accuracy. From there, and 



using purely analytical expressions, we can directly ob- 
tain the second time derivatives of the mass quadrupole 
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and octopole moments and the first time derivative of the 
current quadrupoie moment. Thus, we just need to com- 
pute up to the six additional time derivative for the mass 
moments (i.e., starting from their second time derivative) 
and up to the five additional time derivative for the cur- 
rent quadrupoie moment (i.e., starting from its first time 
derivative). 

Computing numerical derivatives, in contrast to nu- 
merical integration, is a subtle task (see, e.g. [108J). For 
instance, if we consider finite difference formulas for the 
different derivatives, the computation requires the partic- 
ular combinations of the function we want to differentiate 
at points close to the evaluation point, and these combi- 
nations are divided by a power of the offset between the 
different evaluation points. If one chooses the offset to be 
too small, high-order cancellations in the combinations of 
the function evaluations can occur beyond machine pre- 
cision, yielding a meaningless final result. If instead one 
chooses the offset to be too big, the function might be 
evaluated at points where its behavior is very different 
from the one near the evaluation points, which in turn 
can also lead to large errors in the numerical derivatives. 
In many situations one can find an interval of offset val- 
ues in which the high-order derivatives are sufficiently 
accurate, but such an interval depends on the orbit char- 
acteristics and it is not easy to predict. Although we tried 
many different finite difference rules (from rules involving 
a few points to rules involving more than 20 evaluation 
points) , as well as other generic numerical differentiation 
techniques (such as numerical interpolation or Cheby- 
shev differentiation) , a large amount of fine-tuning that 
was difficult to predict seemed essential in all cases. 

The key point to improve the differentiation algorithm 
is to realize that the methods we have just discussed 
are more general than needed for the EMRI problem. 
In the latter, one is always dealing with functions with 
certain properties that can be exploited to construct a 
better numerical differentiation method. The key fea- 
ture is that multipole moments are functionals of the 
trajectories, which are piecewise timelike bounded Kerr 
geodesies, and in turn can be characterized by three fun- 
damental frequencies (in the generic case, see the discus- 
sion in Sec. Ill A). Following [86], we know that a general 



there are also two independent frequencies, but they are 
Q e and fL. 

Our procedure to estimate time derivatives is then to 



functional of Kerr orbits, let us call it f[ip, x^ Wj can be 
expanded in a multiple Fourier series of these frequencies, 
that is 

where (fc, m, n) are integers running from — oo to +oo and 
fkmn are complex coefficients such that f_ k _ m _ n = 
fkmn- There are three special cases in which this ex- 
pansion is simplified: (i) Circular-equatorial orbits; (ii) 
Equatorial noncircular orbits; (iii) Circular nonequato- 
rial orbits. In case (i), the Fourier series contains only a 
single frequency, the azimuthal one fi^. In case (ii), there 
are two independent frequencies, fl r and fL. In case (iii), 



first fit an expansion like that of Eq. ( 97 1 to the multipole 



moments that are required using a standard least-square 
fitting algorithm, and then to estimate time derivatives 
via 

f^,X,m)= E fk,m,ne- i(knr+mQs+nnJt , (98) 

k,m,n 

where 

fk, m , n = H) N (kn r + m n e + nn^f f k ^ n . (99) 

The Fourier fits to the multipole moments, therefore, 
play a crucial role in the accuracy of the high-order 
derivatives. We carry these fits out by evaluating the 
function to be fitted on a certain number of points along 
a geodesic piece of the orbit. As we have already men- 
tion, we can compute analytically the first time deriva- 
tives of the multipoles, so the time-dependent functions 



that we actually fit are: M^(t), M\f h {t), and S^'(t) for 



,(i) 



the radiation-reaction potentials, and also M ijkl (t) and 



ijkl \ 

Sijk(i) f° r the waveforms. The parameters that we need 
to choose for the least-squares fit are: (i) the size in time 
of the interval where we fit the function; (ii) the num- 
ber of points in this interval where the function is going 
to be evaluated (i.e. the number of points to be fitted); 
and (iii) the number of harmonics/frequencies that we 
include in the finite expansion of Eq. ([97]). For choice 
(i), we take a fixed fraction of the shortest orbital period 
(i.e. of the minimum of T r = 27r/f2 r , T e — 27r/J7 e , and 
= 2n/Q^; see Appendix|F). For choice (ii), we use be- 
tween 50 — 500 points, depending on the case we are deal- 
ing with (generic or very particular) and the precision we 
want to achieve. For choice (iii), we use between 2 and 5 
harmonics of the fundamental frequencies. Adding more 
harmonics would increase the accuracy of the derivatives, 
but we have empirically found that 5 harmonics is usually 
sufficient. 

For the practical implementation of the least-squares 
fit we use the GNU Scientific Library |105j . We have 
performed a number of experiments for different types of 
functions (and also for multipole moments of the orbital 
trajectory) and we have found that this technique is very 
robust and provides very high accuracy even for the high- 
est derivatives. For instance, for the sixth time derivative 
we find typical accuracies of one part in 10 5 . Taking into 
account that the magnitude of the time derivatives of the 
multipole moments decreases significantly with the order 
of the derivative, this accuracy is more than enough for 
our purposes. Another important feature of this tech- 
nique is that is has appeared to be quite robust with 
respect to the three choices of parameters we have dis- 
cussed above. 
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V. NUMERICAL RESULTS 

In this section, and in order to illustrate the new kludge 
scheme, we present some numerical results from a nu- 
merical code that we have developed to implement new 
kludges, as well as some comparison with other results 
in the literature. The examples of new kludge evolutions 
shown here are all for prograde orbits but there are no 
obstacle to produce similar results for retrograde orbits. 

Let us first consider our proper use of harmonic coordi- 
nates, an ingredient of new kludges that is very different 
from traditional kludge implementations. In the latter, 
(see, e.g. |39| ) harmonic coordinates are approximated 
via (Euclidean) Cartesian Boyer-Lindquist coordinates, 
i.e. (x BL , y BL , z BL ) = (rsin0cos<^,rsin$sin0, rcos#). 
Such coordinates will differ from true harmonic coordi- 
nates greatly in the strong-field regime. In turn, this will 
modify the resulting trajectories and waveforms, as the 
proper choice of coordinates is crucial in the calculation 
of multipole moments, both in radiation-reaction compu- 
tations and in waveform production. 

Let us then compare how much error is introduced by 
the use of the wrong coordinate system in EMRI wave- 
form construction. For this, we employ circular equa- 
torial orbits as they are simpler when comparing the 
waveform phase. Fig. [T] shows a circular-equatorial in- 
spiral orbit for a system with: a/M # = 0.1, q = 1/10 
and p Q = 10, where p Q is the initial value of the semi- 
latus rectum. That is, both orbits have the same initial 
scmilatus rectum, and hence the same initial energy and 
orbital frequency Q, i. We have chosen such a large mass 
ratio in this case, so that one could see the trajectory 
tracks in the figure. The orbit on the left panel has been 
built and represented using the harmonic coordinates 
whereas the orbit on the right has been built and repre- 
sented using the Cartesian Boyer-Lindquist coordinates 
x BL (these coordinates are used both for the integration 
of the geodesic equations of motion and for the estima- 
tion of the self-force). 

The trajectories in this figure are stopped at the radius 
of the last stable circular equatorial orbit (LSO), given 
by [84]: 



r LSO = M. 1 3 + Z 2 T yj(3 - Z x ) (3 + Z t + 2Z 2 ) j , 

(100) 

where the upper (lower) sign is for prograde (retrograde) 
orbits, and 



Z 1 = 1 + 1 - 
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One must be careful when comparing trajectory tracks, 
since these are completely coordinate dependent (by def- 
inition). Instead, we can compare the number of GW cy- 
cles, which are directly observable. The number of orbital 
cycles in the case of harmonic coordinates is essentially 
double that of the Boyer-Lindquist Cartesian coordinate 
case, and consequently, the same is true for the associated 
waveforms (not shown here). This extremely large dis- 
tance is perhaps a bit of an overestimate, since we consid- 
ered very strong-field EMRIs and used the two different 
systems of coordinates in both the waveform generation 
and the calculation of the multipolar, post-Minkowskian 
self-force. This last point is important because in many 
kludge schemes the coordinates only enter in the tra- 
jectory and in the waveform construction, whereas the 
radiation-reaction part is based on PN results or on BH 
perturbation theory. In any case, this example shows 
that the proper and consistent choice of coordinates can 
play a major role in the final waveforms produced. 

Let us now consider how other choices in the new 
kludge scheme modify the final waveforms produced. 
One approximation one can make to simplify new kludges 
is to ignore the local potentials K^ u and introduced 
in Eqs. pSt-pS]), i.e. to pick K^ v = = 0, or in the 
case that we use the approximate harmonic coordinates 
of Appendix [C] we can use the expansions of Appendix [D| 
for these potentials to any order. Obviously this can 
make a big difference and in our simulations we have al- 
ways used the Kerr local potentials in exact harmonic 
coordinates [Eqs. (54)- (56)]. 

Another approximation one can make is to use just 
the radiation-reaction potential of Burke and Thorne 
[Eq. ( [47| ] (with = 0), instead of the full potentials 
of Eqs. @ and @ (see [HJ [§3). In order to illus- 
trate the waveform difference in this case we have stud- 
ied the inspiral of a system with M = 4.5 x 10 6 M Q , 
a/M, = 0.98 and q = 10~ 5 . In Fig. [5] we show the evo- 
lution of the semilatus rectum for two cases, one with 
p(t = t Q ) = p a = 3 (the plot on the left), being t a the 
evolution initial time, and the other one with p Q — 8. As 
we can see from the figure, the higher-derivative correc- 
tions to the Burke-Thorne radiation-reaction potential 
[Eqs. (44) and (45)] increase the radiation-reaction ef- 
fects, in the sense that p and the other orbital elements 
change more rapidly when these corrections are included. 
Moreover, these corrections are more significant in the 
strong-field region, near the last stable orbit, and less so 
as the distance to the MBH increases. 

We can assess quantitatively the difference due to in- 
troducing corrections to the Burke-Thorne potential by 
looking at the GWs emitted. To that end, we introduce 
the following definition: 



Z2 ~ YM + Zh 

For the evolutions shown in Fig. [T] the LSO r « 5.67Af # 
and r H sa 4.67Af # , consistent with r H ~ r — Af # . 



ihs. 



(101) 



where A 



h\ is the GW amplitude and 

arctan (h x /h + ) is the accumulated GW phase. 
The GW phase difference induced by the presence of the 
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FIG. 1. Circular and equatorial inspiral of a binary system characterized by a/M, = 0.1 and q = 1/10. The inspiral starts 
at (x,y) — (10,0)M.. The figure on the left shows the inspiral using harmonic coordinates and the figure on the right using 
Cartesian Boyer-Lindquist coordinates. 
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FIG. 2. Evolution of the semilatus rectum in circular and equatorial inspirals of a binary system characterized by: M. = 
4.5 x 10 6 M Q , a/M, = 0.98, and q = 10" 5 . The fi gure on the left shows an inspiral that has started with p = 3 and the figure 
on the right an inspiral that has started with p a — 8. While the one on the left stopped near the last stable circular orbit, the 
one on the right was stopped after 0.2 yr of evolution. 



corrections to the Burke-Thorne potential is shown in 
Fig. [3] for the evolution with p a = 8 that corresponds to 
the right panel of Fig. [2] Observe that the GW phase dif- 
ference increases with time to accumulate up to 4.54 rad 
for a total evolution time of 0.2 yr. This means we can 
expect a dcphasing of more than 3 cycles for a 1 yr evolu- 



tion, i.e. the radiation-reaction corrections to the Burke- 
Thorne potential are of relevance for precise and long 
EMRI evolutions. The situation is even more dramatic 
for the strong-field evolution that starts with p Q = 3 (left 
panel of Fig. [2]) , where the evolution with only the Burke- 
Thorne potentials takes 55.2 days while the one with the 
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full radiation-reaction potentials takes 42.8 days before 
reaching an unstable orbit, which translates in a differ- 
ence of 2862.5 GW cycles. 




Fig. [2] we find no difference in the GW phase, but a 3.3% 
(for the case with p Q = 8) and 4.8% (for the case with 
p a = 3) difference in the averaged GW amplitude. In 
this comparison we have used the full radiation-reaction 
potentials for all evolutions. 

Let us now present some results for orbits more generic 
than circular equatorial, since our kludges can easily han- 
dle these as well. Consider first circular nonequatorial or- 
bits, since these are the simplest kind of orbits that allow 
us to study the evolution of the orbital inclination, either 
as described by 6 inc [see Eq. (El )] or by 1 [see Eq. (25 1]. In 



this sense, it is interesting to test the conclusions of [42] 
(see also [109] ) that inclination remains almost constant 
for most such EMRIs and also to compare some quanti- 
tative results in order to assess the present new kludge 
implementation. 

The approximation dijdt = has been considered by 
different authors in order to obtain an evolution equation 
for the Carter constant: 



di 
dt 



= 



C- 2C L 



(102) 
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FIG. 3. Evolution of the GW phase difference between two 
inspirals characterized by: M. = 4.5 x 10 6 M Q , a/M. = 0.98, 
q = 10 -5 and p a — 8. The GW phase difference is computed 



using the formula: A<I? GW (i) = $ 



urke- 
GW 



3 (*)-*c 



where $^ ke " 

uses only the Burke-Thorne potential [Eq. ( 47 )] and $ 



is the GW phase for an evolution that 



p.Full RR 



is the GW phase for an evolution that uses the full radiation- 
reaction potentials [Eqs. |54|-(l56|]. 



Finally, let us look at the choices associated with the 
waveform construction model. In this part of the new 
kludge scheme, we can choose the multipolar order of the 
expansion of the gravitational radiation field as described 



Given that it is not simple to estimate the evolu- 
tion of C within the framework of perturbation theory 
(see [58 l IllOj ). this approximation provides a clear path 
to the construction of generic EMRI waveforms. The PN 
leading-order prediction for the evolution of the inclina- 
tion given by Ryan |111) predicts a grow in the incli- 
nation. Although this prediction is known to overesti- 
mate the inclination growth (see [109] for a discussion), 
it agrees with BH perturbation computations valid in the 
strong field. 

Table [TT] presents some results of the evolution of the 
radius r a and Carter constant C for circular nonequato- 
rial orbits characterized by: 

(i) a/M. = 0.05, rjM. = 100, and 1 = 60.0 deg; 



in Sec. (HIE} Three possibilities present themselves here: ( U ) a l M ' = °' 95 ' r °/ M ' = 100 < and L = 60 - 05 de & 
(i) Quadrupolar waveforms; (ii) Octopole-Quadrupolar (m) = Q ^ /M> = ^ and L = 6Q U deg; 

waveforms; and (iii) Hexadecapole- Octopolar waveforms. 

Case (i) is the lowest-order approximation and consists of (i v ) alM % — 0.95, r /M. = 7, and 1 = 60.43 deg . 
considering only the mass quadrupole term in the wave- 
form expansion. This is equivalent to considering up to 
second time derivatives of the trajectory, i.e. up to accel- 
erations. Case (ii) accounts for the next order multipole, 
that is, the mass octopole and the current quadrupole. 
This is equivalent to considering up to third time deriva- 
tives of the trajectory, i.e. up to the jerk: f = d 3 x i /dt 3 . 
Case (iii) adds one more multipole, that is, the mass hex- 
adecapole and the current octopole. This is equivalent to 
considering up to the fourth time derivatives of the tra- 
jectory, i.e. up to the snap: s l = d 4 x l /dt 4 . 

In general, adding more multipoles does not affect the 
GW phase but it introduces amplitude corrections that 
depend on the inspiral character. For instance, compar- 
ing the quadrupole waveforms with the hexadecapole- 
octopole waveforms for the evolutions corresponding to 



These evolutions use Eq. (|36j) and the formulas of Ap- 
pendix [G] Observe from the table that the new kludge 
results are in very good agreement with Teukolsky like 
evolutions, even though the former has not really been 
refined or optimized. 

The local nature of new kludges can also be appreci- 
ated in the local, nonuniform temporal changes of the 
inclination angle, i.e. di/dt is not constant in time but 
oscillates with the orbital period T e (see Table 
value of the fundamental frequencies f2 
cases of Table |n| . In contrast, we find that dr /dt is ap- 
proximately constant within orbital time scales. We illus- 
trate these facts in Fig. [4] where we compare the evolution 
of both quantities for a total time of At — 800 M. and for 
the third case of Table [TH characterized by a/M # = 0.05, 
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a/M. rJM. i (deg) 


Quantity 


PN Teukolsky New Kludge New Kludge 
(Ref. [55]) (Ref. [43;) Burke-Thorne Full RR 


0.05 100 60.00 


q- 1 {drjdt) 
{M./q) (dt/dt) 


-1.2797 x 10" b -1.2676 x 10" 5 -1.1700 x 10~ b -1.1708 x 10" 5 
7.0439 x 10" 12 6.6936 x 10~ 12 6.1597 x 10~ 12 6.5089 x 10~ 12 


0.95 100 60.05 


q- L {drjdt) 
{M./q) {dt/dt) 


-1.2733 x 10" 5 -1.2610 x 10" 5 -1.1622 x 10"~ b -1.1634 x 10 _b 
1.3389 x 10" 10 1.2040 x 10~ 10 1.1628 x 10~ 10 1.2273 x 10" 10 


0.05 7 60.17 


q- 1 {drjdt) 
{M./q) (dt/dt) 


-3.6762 x 10~ 2 -1.0964 x 10 _i -8.3338 x 10^ -8.9633 x 10~ 2 
1.5867 x 10" 5 1.0875 x 10" 5 9.6088 x 10" 6 1.6233 x 10~ 5 


0.95 7 60.43 


q- L (drjdt) 
{M./q) (dt/dt) 


-2.7499 x 10~ 2 -4.6574 x 10~ 2 -3.4547 x 10~ 2 -3.6825 x 10~ 2 
3.0806 x 10~ 4 1.2073 x 10" 4 1.6023 x 10" 4 2.5962 x 10" 4 



TABLE II. Evolution of the radius, r a , and inclination angle, t, of circular nonequatorial orbits characterized by (columns 1st 
to 3rd from top to bottom): (i) a/M. = 0.05, rJM. = 100, and i = 60.0 deg; (ii) a/M. = 0.95, rJM. = 100, and l = 60.05 
deg; (iii) a/M. = 0.05, rJM. = 7, and t = 60.17 deg; (iv) a/M. = 0.95, rJM. = 7, and t = 60.43 deg. Column 5th gives the 
value of q~ (drjdt) and (M./q) (dt/dt) obtained from the PN calculations of Ryan |89| ; column 6th gives the value obtained 
by Hughes ,43, solving the Teukolsky equatio n; an d columns 7th and 8th are t he v alues obtained with our current new kludge 
implementation using the Burke-Thorne [Eq. (47 1] and the full [Eqs. (441 and (45 1] radiation-reaction potentials respectively. 



a/M. 


r /M. 


l (deg) 






M.Q^ 




0.05 


100 


60.00 


9.9992 x 10" 


-4 


1.0000 x 10" 


-3 


0.95 


100 


60.05 


9.9856 x 10" 


-4 


1.0004 x 10" 


-3 


0.05 


7 


60.17 


5.3776 x 10" 


-2 


5.4065 x 10" 


-2 


0.95 


7 


60.43 


4.9813 x 10" 


-2 


5.4537 x 10" 


-2 



TABLE III. Values of the fundamental frequencies Q g and 
(see Appendix [F] for details) for the four cases shown in 
Table [TT] The values of the frequencies are normalized with 
respect to the MBH mass M. . 



rJM, = 7, and i = 60.17 deg. The results of Table [TT] 
can be obtained via a simple linear regression of the evo- 
lution of i over a number of orbital periods (we quote the 
slope as the value of dt/dt in Table |TT|) . 

Let us comment further on the results of Table [TT| clas- 
sifying them into weak-field calculations [cases (i) and 
(ii)] and strong-field ones [cases (iii) and (iv)]. First, as 
it was already discussed, the new kludge results in the 
weak field do not differ much using the Burke-Thorne 
or the full radiation-reaction potentials, whereas in the 
strong field the differences are significative. Second, con- 
sidering the absolute value of the numbers quoted in the 
table, new kludges tend to always underestimate the rate 
of decay of the radius of the circular orbit, r a , with re- 
spect to the Teukolsky results of [23] ■ This is in contrast 
with the PN results of [55], which overestimate dr a /dt in 
the weak field and underestimate it in the strong field. 

Regarding the evolution of the inclination angle t the 
situation is a bit different, although we can still see sim- 
ilar differences between using the Burke-Thorne and the 
full radiation-reaction potentials in the strong-field com- 
putations. For the weak-field computations, the differ- 
ences are significantly bigger for the evolution of t than 
for the evolution of r a . To put these numbers in context, 
the errors in the new kludge estimations of the rates of 
r Q and t, as compared with the Teukolsky estimations, 
are: for dr a /dt, 7.63 — 7.83% for the weak-field com- 
putations and 18.24 — 25.82% for the strong-field ones; 




200 400 600 



Time (M.) 

FIG. 4. Evolution of the radius r (upper plot) and inclination 
t (lower plot) for the circular nonequatorial orbit character- 
ized by a/M. = 0.05, rJM. = 7, and t = 60.17 deg. This is 
the third case in Table [TT1 The total time of the evolution in 
units of the MBH mass is 800M. . 



for dt/dt, 1.93 — 7.97% for the weak-field computations 
and 11.64 — 115.04% for the strong-field ones. In the 
weak field, the new kludge computations, again assum- 
ing that the Teukolsky computations are the correct ones, 
in general do worse than the PN ones for dr Q /dt but they 
do better for dt/dt, whereas in the strong field the new 
kludge computations seem to be better than the PN ones 
in general. Of course, this assumes that the Teukolsky 
waveforms are exactly correct, which is not the case ei- 
ther. A more detailed comparison between Teukolsky 
and new kludge results will be carried out elsewhere. 
Let us finish by considering an example of generic, 
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eccentric and inclined orbital evolutions with the new 
kludge scheme. For such orbits, all the constants of mo- 
tion/orbital elements change in time, plotted in Fig.Jsjfor 
a sufficiently long time, including many orbital periods. 
The main figure hides the changes in the orbital time 
scales, so we have included subplots where we can see 
that the evolution details scale with the orbital periods. 
As we have already discussed above, this is a consequence 
of using a local in time self-force. These effects are not 
present in evolutionary schemes based on flux averaging 
over a certain number of orbital periods (essentially all 
other models currently used). 

As we can see in Fig. [5j the local evolution of the dif- 
ferent quantities presents slightly different patterns. The 
difference in these patterns are enhanced if one chooses 
a more extreme orbit (more eccentric and more strong 
field). We can also appreciate that looking at the global 
evolution (over the long time scale) all quantities decay 
in time except for the inclination, which grows in time. 
However, if one looks at the details of the evolution over 
the orbital periods, one can see (eg. for the eccentric- 
ity) that it can locally grow in time although the global 
tendency is to decay. Therefore, the new kludge scheme 
leads to quite rich evolutionary patterns due to its local 
in time character, which also makes it a very valuable tool 
to investigate questions like the appearance of transient 
resonances in EMRIs [73] . 

To illustrate the multipolar waveforms produced by the 
new kludge scheme for general orbits, Fig. [6] shows short 
fragments (for the sake of clarity) of the GWs associated 
with the last evolution of Fig. [5j i.e. for the inspiral of a 
system with M. = fO 6 M , a/M, = 0.98, and q = f(T 5 , 
and with initial orbital elements (p, e,i) = (7,0.6,57.39 
deg). The waveforms polarizations shown in Fig.|6j which 
correspond to an observer along the spin axis, present the 
richness of EMRI GWs for eccentric and inclined orbits. 



VI. CONCLUSIONS AND DISCUSSION 

We have introduced the new kludge scheme to model 
the dynamics and the GW emission of EMRIs, which 
in principle could also be used for intermediate-mass 
ratio systems. This scheme combines ingredients from 
the multipolar, post-Minkowskian formalism and black 
hole perturbation theory to evolve a nongeodesic word 
line (with respect to the geometry of the binary's large 
component, assumed here to be a MBH) and construct 
waveforms. The orbits are built as a sequence of local 
geodesies whose orbital elements evolve according to a 
local self-force that we approximate via a multipolar, 
post-Minkowskian expansion. The leading order term 
of this self-force corresponds to the well-known Burkc- 
Thornc radiation-reaction potential. We have seen that 
a crucial ingredient in this construction is the mapping 
from Boyer-Lindquist coordinates, in which we integrate 
locally in time the orbits to the harmonic coordinates 



required both by the multipolar, post-Minkowskian self- 
force and for the GW multipolar expansion. Once we 
have trajectories in harmonic coordinates, it is straight- 
forward to build waveforms. 

In practice, the implementation of the new kludge 
schemes requires a number of numerical/analytical in- 
gredients. First, one must numerically integrate (locally 
in time) geodesic equations, for which we use appropriate 
angle variables that avoid turning points. From the inte- 
gration of these equations, we can use analytical formu- 
las to map the trajectories (and associated velocities and 
accelerations) to harmonic coordinates. Then, informa- 
tion from small fragments of the geodesic orbit (smaller 
than the fundamental orbital periods) is used to build 
the different multipole moments and compute their time 
derivatives. We have found this to be the most chal- 
lenging point of the implementation as we need up to 
eight-order time derivatives of the multipoles, in partic- 
ular of the mass quadrupole. Numerical differentiation 
is much more complicated than numerical integration, 
since the choice of the offset in finite difference formu- 
las is crucial to obtain correct results. By experimenting 
with different numerical techniques we have concluded 
that general differentiation techniques (i.e. valid for any 
differentiable function) are quite difficult to implement 
successfully (specially for the high-order derivatives re- 
quired), since the derivatives require precise fine-tuning 
that depends on the fundamental frequencies of motion. 
The way out has been to use a numerical method that 
takes into account what we know analytically about lo- 
cal motion: the fundamental frequencies. More specif- 
ically, we find that fitting a truncated multiperiod (the 
fundamental periods, which we also obtain numerically 
in terms of elliptic integrals) Fourier series to fragments 
of the local evolution using standard least-squares meth- 
ods provides the accuracy that we require for a precise 
estimation of the radiation-reaction potentials and the 
self-force. With this, we can then evolve the constants 
of motion, thus mapping from the initial geodesic to new 
ones. Such a transition requires the mapping of the new 
constants of motion (E, L Z ,C/Q) to the new orbital el- 
ements (e,p, i/0 im> ), which can be done using analytic 
formulas. Repeating this procedure iteratively in time, 
we build the inspiral trajectory or worldline, from which 
we then finally construct the waveforms. 

Apart from the exact harmonic coordinates of Ding 
(see Ref. |98j ) that we use in our implementation, we have 
also provided approximate harmonic coordinates based 
on the construction of ACMC coordinates (see Appen- 
dices [C] and [D]) . The construction of these coordinates is 
a useful exercise that can be used in scenarios where the 
big component of the binary is not described by the Kerr 
geometry but by something else either related with other 
theories of gravity or with the idea that MBH at galactic 
centers may be exotic matter configurations rather than 
black holes. 

The results presented here serve as a first introduction 
and proof-of-principle of the new kludge scheme. Much 
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FIG. 5. Evolution (for a total time of 10 -2 yrs) of an eccentric and inclined inspiral of a system characterized by: M. = 
10 6 Mq, a/M. = 0.98, and q = 10 -5 . The plots show the evolution of the following quantities: Energy E (top left), actually 
10 4 [-E — E(t a )], where E(t a ) — 0.9575513; angular momentum along the spin axis, L z (top center); Carter constant C (top 
right); semilatus rectum, p (bottom left), with p(t a ) = 7; eccentricity, e (bottom center), with e(t ) = 0.6; and inclination angle 
l (bottom right), actually 10 4 [t — i(t )], where t(i ) = 57.39 deg. All plots contain subplots where the detailed evolution during 
a few orbital periods is shown. 




2000 4000 6000 

Time (sec) 

FIG. 6. Fragments of the GWs emitted from an eccentric and 
inclined system [with initial orbital elements (p,e,i)(t ) = 
(7,0.6,57.39 deg)] characterized by: M. = 1O 6 M , a/M. = 
0.98, and q — 10 . The evolution of the constants of motion 
and orbital elements of this system are shown in Fig. [5] The 
solid line represents the + polarization, h + (t), and the dot- 
ted line represents the x polarization, h x (t), as seen by an 
observer located along the spin axis. 



work remains to be done to improve and to validate the 
method in order to obtain GWs to the level of accuracy 
required for LISA data analysis. One clear way of im- 
proving the scheme is to use better expressions for the 



different multipole moments that go beyond the present 
leading-order approximation. One caveat of doing this 
is that the PN corrections required are in general not 
known for generic spinning binaries. In any case, the 
currently known PN corrections will certainly improve 
the accuracy of the new kludge evolutions. Another way 
to improve the scheme would be to introduce conserva- 
tive corrections to the background, for instance as it is 
done in the EOB formalism [48 . 

A more detailed and exhaustive validation of the new 
kludge scheme would include, as a first step, a compari- 
son of the evolution of the constants of motion with the 
fluxes associated with them that are calculated by solv- 
ing the Teukolsky equation [42j 03] . As we have already 
mentioned above, one has to be careful in doing this com- 
parison given that the latter employs averages over sev- 
eral cycles, while the new kludge fluxes are computed 
locally at the SCO's location. Once the fluxes have been 
validated, one should compare the waveforms themselves. 
An overlap study would determine the level of agreement 
between them. 

Another interesting aspect that we can test is whether 
the new kludge scheme can be used for IMRIs and even 
for systems with moderate mass ratios. Recently, com- 
parisons between self-force, PN, and numerical relativity 
computations have arrived to the conclusion that by re- 
placing the mass ratio q by the symmetric mass ratio 77, 
the self-force predictions compare quite well with numeri- 
cal relativity and PN predictions in the comparable-mass 
range (331 135] . These results may allow for a simpler de- 
scription of IMRIs, which otherwise would require either 
long full numerical computations or higher-order pertur- 
bative computations, or a combination of both. In the 
new kludge scheme, q — > 7/ is actually mostly built in 
already, since the mass ratio information enters through 
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the definition of the multipole moments [see Eqs. (48 1, 
and (86)]. For these, we already use gen- 



(49 



era 



(85), 



binary expressions, instead of effective one-body ones 
with mass qM 9 . Therefore, it would be interesting to see 
whether future improvements of the new kludge scheme 
can produce reasonably good results for systems other 
than EMRIs. 

The ability to compute approximate local quantities 
at the location of the SCO, in particular the multipo- 
lar, post-Minkowskian self-force, makes the new kludge 
scheme an interesting tool to study certain local behavior 
believed to exist in EMRIs. For example, Flanagan and 
Hinderer [73j have recently reported that certain rapid 
changes in orbital elements can arise for generic EM- 
RIs when the orbital frequencies become commensurate. 
During these rapid changes, it has been postulated that 
the EMRI waveforms might suffer a "glitch," not unlike 
those observed in pulsar astronomy. Questions remain 
though as to the exact nature of this effect. One could 
compare the effect of these glitches in waveforms as com- 
puted from an averaged-scheme (such as the Teukolsky 
one) and a local one (such as the new kludge approach). 
These, and other studies, would provide one more piece 
to the EMRI puzzle. 
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Appendix A: Useful formulas for the Computation 
of the Self- Acceleration 



In this appendix we define the quantities relevant to 
the multipolar, post-Minkowskian radiation-reaction ac- 



celeration of Eq. (|57|), where the pieces Aa and Aa' are 



I (2) 



given in Eqs. (58) and (60) respectively. We also need 



the expression of the projector orthogonal to the SCO 
four-velocity. This quantity can be written in terms of 
the potentials and as follows: 



1"P 



+ TW 3 . 



p. 



(Al) 
(A2) 



where the T factor was already given in Eq. ( 28 ) , but it 



can also be rewritten, in terms of the K„ v potentials, as 

1 



r 



i 



K - 2K,v l 



(A3) 



where we can see how the relativistic T factor is modified, 
with respect to the Special Relativity form T = (1 — 
u 2 ) -1 / 2 , due to Kerr local potentials. 

The first piece of the self-acceleration, A a , is deter- 



JRR 

' fiva 



[Eq. ( 59 1] . The time and spatial 



mined by the object Q 
components of Aa [see (61 )] are given by 

-4 RR = (1 - v 2 ) d t V RR + 2v%V ia . 



(A4) 



1 



diV RR + 2 Vi d t V RR + 2v i v>d j V RR 

- 4au RR • 



(A5) 



The second piece of the self-acceleration, A a ' , is lin- 
ear in the radiation-reaction potentials V RR and V RR [see 
Eq. ( 60 )] . The coefficients of the linear combinations of 
these potentials in the time and spatial components of 

(2) \ V \ V 

A a [Eqs. (62) and (63) respectively] have the following 
form: 



B RR 

^RR 



VI, 



= Q l d t K , 

= 2(1- Q) v%K + \Q'Vi),K, , 

= 2(l-Q)v ij d i K j 

- QV' fc (2d u K k)l - diK jk ) , 

= -2(S ij + Q ij )djK, 

= iQ'v^djK + 8 (6 tk + Q lk ) v 3 d {] K k] 

= AQV k 8 3 K k 



2 5 l 



(2d u K, 



k)l 



diK 



(A6) 
(A7) 

(A8) 
(A9) 
(A10) 

(All) 



All contractions in these expressions are to be performed 
with the Kronecker delta, 5^ = 5^ = diag(l,l,l), and 
we recall that Latin indices in the middle of the alphabet 
stand for spatial coordinates only (thus, all tensors are 
purely spatial). 

For the purposes of the new kludge scheme we need to 
compute these local potentials associated with the Kerr 
metric (K H , Kf, and Kfj) and its inverse (Q H , Q l H , and 
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Qh) in harmonic coordinates. They are need for the or- 



and (A2)], for the V factor 



thogonal projector [Eqs. (Al) 
[Eq. (28)], and for the coefficients of the second piece 

fn\ 

of the acceleration, A a ' [Eqs. (A6)-(All)]. Moreover, 
for these coefficients we also need to compute the first 
spatial derivatives of the metric potentials (d k K H , d k Kf, 
and d k Kfj). This can be done in an efficient way by using 
the following relations: 



Q K _ K,BL BLpM ^BL 



d k Kf 



f) T n a rr, 

^K,BL BLpM BL L/U, B: 

=M(t m)n Q x k Q x i 



(A12) 



three-dimensional flat-space vector algebra: a position 
vector r H = (a; H , ?/ H , z H ), whose spatial norm (under the 
flat three-dimensional metric) was defined in Eq. (74), 
and the reduced spin vector s which, for consistency with 
the choices made in the paper we take it to be aligned 
with the z-axis: s = (0,0, a). Nevertheless, one does not 
need to assume this specific form for s, as the expressions 
we present in this Appendix are invariant under the ac- 
tion of the rotation group. We can then use the above 
expressions for arbitrary directions of the spin angular 
momentum of the BH. Note from Eq. (B2| that we al- 



ways have 



> 



ft > < 



and clearly the equality 



holds for the Schwarzschild case (a = 0). These inequal- 
ities translate in the following inequalities for r: 
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2 &tr. 
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The advantage of these expressions is that they only in- 
volve the Jacobian and Hessian of the coordinate trans- 
formation between Boyer-Lindquist and harmonic coor- 
dinates, and the Christoffel symbols in Boyer-Lindquist 
coordinates. The latter are relatively simple functions 
and have to be computed anyway for other purposes 
(e.g. to compute the accelerations of Boyer-Lindquist co- 
ordinates) . 



Appendix B: Harmonic Coordinates for Rotating 
Black Holes 



In this appendix we describe more properties of the 
set of harmonic coordinates for rotating BHs employed in 
this paper. Eqs. (|67|-([70|) and Eqs. ([7l)-(|73| provide the 



explicit expressions for the transformation from Boyer- 
Lindquist to harmonic coordinates and the inverse trans- 
formation. From these transformations, we find that the 
relation between the Boyer-Lindquist radial coordinate r 
and the spatial harmonic coordinates (x H ,y H ,z H ) is: 



vl 



(r-M.) 2 + , 



(r - M.) 



1 



(Bf) 



which resembles the relation between r and the spatial 
Kerr-Schild Cartesian coordinates (see, e.g. |112j ). This 
relation allows us to write r as a function of r H as given 
in Eq. (72), which we find convenient to rewrite as r = 
M, + g, where 



Q : 



r\ - a? + y/(rl - a 2 ) 2 + 4(s • r H f 



(B2) 



and where the dot here denotes the flat-space three- 
dimensional scalar product. Moreover, we have here in- 
troduced the following notation in the spirit of the usual 



M. > r > 



(B3) 



An important ingredient of the construction of the har- 
monic coordinates is the angle function $(r), which was 
already given in Eq. ( 75 1 . However, its proper definition 
is actually [98 



$(r) = - / 

J r 



dr' 



aMl 



A(r')[A(r') + Af 2 ] ' 



(B4) 



where we recall that A(r) = r 2 — 2M,r + a 2 . An alterna- 
tive expression for <f>(r) can be obtained by simplifying 
Eq. |75l) into 



* / x t fr-M. 

<J>(r) = arctan 

2 \ a 



: In 
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Two important relations to evaluate the coordinate 



transformation of Eqs. (67)-([70j) are: 



COS $ = 



1 + O2 (r)] { 1+( r^) 2 }' 



(B6) 



sin $ = 



r-M. 



fi(r) 



^[l + 2 (r)]{l + (^) 2 } 



(B7) 



where £l(r) is given in Eq. (76). The explicit form of the 



components of the Kerr metric and its inverse in har- 
monic coordinates is given in Appendix A of |98j . These 
expressions are quite lengthy and, in that form, not very 
efficient for numerical computations. In what follows, we 
present an alternative method to compute the compo- 
nents of the transformed metric in an efficient way. At 
the same time, we introduce expressions that are invari- 
ant under the rotation group, in the sense that they do 
no longer assume the BH spin axis is aligned with the 
z„-axis. 
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Jacobian 



components of the metric in harmonic coordinates are 



We start by presenting the expressions of the com- 
ponents of the Jacobian associated with the coordi- 
nate transformation between Boyer-Lindquist and har- 
monic coordinates. In other words, we compute the 
matrices J| L = D(t,r,9, <j))/ D(t a , x li ,y H , z H ) and Jg L = 
D(t H ,x H ,y H ,z H )/D(t,r,6,(j>), which obviously are in- 
verses of each other: J| L • J| L = J BL • J" L = Id 3 , where 
Id 3 is the identity matrix in three-dimensions. The com- 
ponents of the Jacobian J| L are 
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where unlabeled coordinates on the right-hand sides 
stand for Boyer-Lindquist coordinates as usual. Simi- 
larly, the contravariant components of the metric in har- 
monic coordinates are 
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where the x denotes the usual flat-space three- 
dimensional vector product. The components of the in- 
verse Jacobian, i.e. J| L , are 
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2. Kerr Metric in Harmonic Coordinates 

The expressions above can be used to compute the 
transformed Kerr metric and its inverse in harmonic co- 
ordinates systematically and efficiently The covariant 



It is easy to see that, under transformation of the three- 
dimensional rotation group, the covariant and contravari- 
ant (i H ,i H ) components transform as scalars, the co- 
variant and contravariant (i H ,r^) components transform 
as vectors, and the covariant and contravariant (r l H7 r^) 
components transform as 2-rank tensors. The expres- 
sions of all these components are invariant. 

The computation of the components of the Kerr met- 
ric and its inverse in harmonic coordinates requires ex- 
pressions for the Boyer-Lindquist components of the Kerr 
metric in harmonic coordinates. The covariant compo- 
nents arc 
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Hessian 



Let us now present explicit expressions for the Hes- 
sian of the transformation between harmonic and Boyer- 
Lindquist coordinates, needed for the computation of the 
accelerations in harmonic coordinates. Also, for the com- 



putation of the derivatives of the metric in harmonic co- 
ordinates we need the Hessians of the Boyer-Lindquist co- 
ordinates with respect to the harmonic ones. More specif- 
ically, for the computation of the coefficients of the sec- 



(2) 

ond piece of the acceleration, AX' [Eqs. (|A6[)-(|A11[)]. We 



give here these expressions using the three-dimensional 
vector notation introduced above. The expressions for 
the Hessian matrix d 2 r^/(dx l BIj dx 3 Blj ) are: 
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and the expressions for its inverse d a;| L /(9r^ dr^) are: 
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where </>r H r H is the following symmetric matrix 

V * * o 



and where denotes symmetric tensor product: aQb 
(a<gib + b<E)a)/2. 
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Appendix C: ACMC and approximate Harmonic 
Coordinate Systems 

In this Appendix, we construct an ACMC-6 coor- 
dinate system and a system of approximate harmonic 
coordinates associated with these ACMC-6 coordinates 
that we compare with the exact harmonic coordinates 
presented in Sec. HID Thorne [57] constructed an ex- 
plicit map between Boycr-Lindquist coordinates (xg L ) = 
(i, r, 9, 4>) and ACMC-2 for a Kerr BH with mass M, 
and Kerr spin parameter a. The extension of this map 
from Boyer-Lindquist coordinates to ACMC-6 coordi- 
nates, (iZacmc) = (T, X, Y, Z), is given below for the first 
time: 
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are three-dimensional 



have introduced the following notation in the Euclidean 
vector calculus style: r AfJ = r l Ali = (x AH ,y A „,z 

r AH = (-^ah + Hah + z ah) 

6 AH = arccos(s-n AH ), y AH 
(spatial), orthogonal vectors in the plane orthogonal to 
the spin axis. 

These approximate harmonic coordinates are not 
only ACMC-6 but also harmonic up to terms of or- 
der 0(M 6 /r 6 H ), as one can check by verifying that 

K,AH 





9rw 
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3K.AH &K,An) 



0, where g a p is the transformed 
Kerr metric ^]and <9* H denotes partial differentiation with 
respect to these coordinates. 

One should note that this coordinates are only pseudo- 
harmonic, in that they do not respect the harmonic co- 
ordinate condition to all orders in M./r AH . Therefore, 
these coordinates differ somewhat from the transforma- 
tion in Eq s. (| 67|)-(70). We can show this by expanding 



Eqs. (71|-(|73[) in M./r H < 1 and a/r H < 1: 
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where we have introduced the shorthands R = (X 2 + 
Y 2 + Z 2 ) 1 ' 2 , and 6 = arccos(Z/.R). 

Harmonic coordinates are guaranteed to be ACMC- 
oo (see [UJ), but the converse is not necessarily true: 
ACMC coordinates of order N are not necessarily har- 
monic. Nonetheless, one can enhance the ACMC-6 co- 



ordinate transformations of Eqs. rtClHC4) to construct 



a map from Boyer-Lindquist coordinates to approximate 
harmonic (AH) coordinates (x Ali ) = (t 
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[1 - 15(5 • n H ) 2 + 35(s • n H ) 4 - 21(1 • n H ) 6 ] , 



2r 2 1 



3a^ 
4r 2 
5a 4 



24r 4 
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j [l-4(s-n H ) 2 ] [3-4(|.« H ) 2 ]j ,(CU) 



aM 2 

0[a 2 -5a 2 (|.n H ) 2 



2Mf 



(C12) 



where we have used a similar vector notation as 



Eqs. (C5|-(C8| with 9 H — arccos(s • n H ) and tf> B = 
arctan(y H • n H /x n ■ n H ). As we can see, disagreements 
arise in the radial transformation at order C(a 4 ), but 
this difference is proportional to a monopole (without 
angular dependence) and a quadrupole term (quadratic 
in angular dependence), which do not modify the ACMC 



arctan 



(C8) 



where s is a unit vector (in the Euclidean sense) along 
the MBH spin axis which, according to the conventions 
of this paper, has components s — (0, 0, 1). Moreover we 



1 This transformation is slightly different from that found in 1 1131 ■ 
for zero integration constants. One can show by direct evalua- 
tion, however, that both transformations lead to harmonic coor- 
dinates that are also ACMC. 
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condition at octopole order. Similar disagreements arise 
also in the angular sector of the transformation. 

We have thus shown that the ACMC coordinate map 
we found in Eqs. (C5)-(C8) reproduces the main ingredi- 
ents of the full harmonic coordinate transformation. In 



fact, since we have shown that Eqs. (C5|-(C8| lead to a 



metric that satisfies the harmonic coordinate condition, 
we can infer that the difference between this equation 



and Eqs. (71)-(73l must amount to a refinement of the 



coordinate system. 



Appendix D: Far-Field Expansion of the Kerr 
Metric in Approximate Harmonic Coordinates 

In this appendix we expand the Kerr metric in the far- 
field using the approximate harmonic coordinates x Ali 
of Appendix [C] [the coordinate transformations from 
Boyer-Lindquist coordinates are given in Eqs. ( C5 )-( C8 )] . 
Then, in transforming the Kerr metric from Boyer- 
Lindquist coordinates to these approximate harmonic co- 
ordinates we assume M m /r AH <C 1 and a/r AH <C 1. We 
give the expressions of the Kerr metric and its inverse 
in terms of the local potentials K a g and Q a P. The far- 
field expansions of the K„ v potentials in the approximate 
harmonic coordinates of Eqs. ( C5 )-( C8 ), K£%, is given by 
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The far-field expansions of the potentials in the ap- 



proximate harmonic coordinates, of Eqs. (|C5|-(|C8|), Qah, 



is given by 
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where we recall that the symbol x refers to the Euclidean 
vector product and the dot product to the Euclidean 
scalar product. Moreover, n AH = x All /r Ali and indices are 
raised and lowered with the flat metric. We have checked 
that the above metric satisfies the differential harmonic 
coordinate condition to the order of approximation. That 
is, we have checked that d Aii (^J— g KjAH g^n) = 0, where 
here g^ H is the expanded Kerr metric in the approxi- 
mate harmonic coordinates as determined by the expan- 
sions of the potentials above, and d AIi denotes partial 
differentiation with respect to the approximate harmonic 
coordinates (x AH ). 



Appendix E: Relations between different orbital 
par ameterizat ions 

Geodesic orbits in Kerr spacetime are fully determined 
by the three constants of motion X A — (E, L Z ,C or Q). 
In this work, we restrict our attention to bounded mo- 
tion (E 2 < 1, see |114j ). Then, it is also very useful 
to characterize the orbit in terms of the orbital param- 
eters O a = (p,e,i or 9 inc ), which provide more trans- 
parent geometrical information about the properties of 
the orbit. Both sets of parameters are important, the set 
I A is more adapted to the separation of the geodesic 
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equations and to the radiation-reaction computations, 
whereas the set O a is better in terms of orbit character- 
ization. Therefore, it is important to know the relation 
between these two sets of constants of motion and how to 
map them. This is specially important when radiation- 
reaction changes these parameters and we need to know 
the new values of O a once the new values of I A have 
been computed using Eqs. p0|-(33 1. The relations we 



present below follow from developments in 

Let us first consider the case in which the set (p, e, inc ) 
is known, i.e. we parametrize the orbit in terms of these 
parameters and we need to find the parameters X A for 
evolving the equations of motion. We mainly use 9 ilic 
instead of i for convenience, but we shall also present 
formulas related to the inclination angle i. The goal is 
to find the mapping from (p, e, 9 inc ) to (E,L Z ,C or Q). 
First of all, given inc , the minimum value that 9 can 
take, min , is [see Eq. 



--sign(Zji 



(El) 



where the sign of L z determines whether the orbit is pro- 
grade (positive) or retrograde (negative) . The sign of L z 
is encoded in inc as follows: If < inc < tt/2, then 
sign(L 2 ) = 1; if —tt/2 < inc < 0, then sign(LJ = -1. 
The particular case inc — is singular in the sense that 
both signs are possible for L z . 

Since 9 = . is a minimum, we have 0(0 mln ) = 0, 



which means that the right-hand side of Eq. ( 20 ) has to 
vanish at min , from which we obtain an expression for C 



C 



1 - z_ 



a 2 {l-E J 



(E2) 



where z_ is given in Eq. ( 92 1 . Now that we have an 



expression for C in terms of E and L z , let us find expres- 
sions for (E, L z ) in terms of p and e. This can be done 
from the analysis of the radial motion, and using the ex- 
pressions of the apocenter and pericenter radii, r apo and 
r , [Eq. p2|], in terms of (p, e). These values of r are 



cxtrcma and hence, the right-hand side of Eq. ( 19 ) van- 
ishes at r = r apo and at r = r pcri , leading to two equations 
for the three unknowns (E, L z , C). The Carter constant, 
however, is given in terms of (E,L Z ) in Eq. (E2|, which 



then leads to two equations for two unknowns, (E,L Z ). 
These equations take the following polynomial structure 



ajE 2 + 2f3jEL z + ll h\ + A 7 = , 



(E3) 



where the subindex / stands for apocenter or pericenter 
and where the coefficients a Iy f3j, 7 j, and A J; in the case 
of noncircular orbits (r wJ ^ ?" apo )j are given by 

cxj = (rj + a 2 ) (r 2 + a 2 z_) + 2A/.r 7 a 2 (l - z_) (E4) 
Pj = -2M.Tja , (E5) 

2M,r 7 ] , (E6) 
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\r 2 T +a 2 z_ 



1 - z_ 

X I = -(rj+a 2 z_)A(r I ) 

= -(rj+ a 2 z_) (rj - 2M.rj + a 2 ) 



(E7) 



In the case of circular orbits these two relations are ex- 
actly the same and we need an extra equation. This 
equation comes from the fact that in the circular case 
dr/dr must always vanish and hence also the radial 
derivative of Eq. ( 19 ) must vanish. Then, the first equa- 



tion in the circular case is given by Eqs. (E4|-(E7| with 
rj = r Q = const., and the coefficients of the second one 
are given by: 

* 2 = 2r D (r 2 a + a 2 



a 2 

72 = 

A, = 



) -a 2 (r a — M.) (l — z_) 



-aM. , 
r -M. 

1 - z_ 



-r A(r )-{r -M,)(r 2 +a 2 z_) 



(E8) 
(E9) 

(E10) 
(Ell) 



One can then combine the two expressions in Eq. ( E3 ) 
to eliminate one of the two unknowns. For instance, we 
can eliminate L z and obtain an equation for E, which 
can be written in the following form: 

([a, 7 ] 2 + 4[a,/3][ 7 ,/?])£ 4 

+ 2 ([a, 7 ] [A, 7 ] + 2[ 7 , ft [A, /?]) E 2 + [A, 7 ] 2 - 0(4312) 

where the notation [*,*] denotes the following antisym- 
metric product: 

[11,121 = n n . -n n , (Ei3) 

L i J apo pen pen apo ' V / 

and the subscripts denote at which radii we have to eval- 
uate the quantity (e.g., I4 apo = II(r )). Eq. (E12| is 



bi-quadratic in E, which means that there are two solu- 
tions for E 2 and from each of these, there are two values 
of E, one positive and one negative, that are related by 
time-inversion. From the two solutions for E 2 , the larger 
one corresponds to retrograde orbits, while the smaller 
one corresponds to prograde ones. Given E 2 , we can 
then find L z through 

L 2 = ^([a,/3]S 2 + [A,/3]) , (E14) 

where the positive solution correspond to prograde or- 
bits and the negative one for retrograde orbits. Fi- 
nally, the Carter constant C is given by Eq. ( E2 ) and 



Q from Eq. (17 1. Once we know the constants of motion 
(E,L Z ,C) we can also find the inclination angle i from 



Eq. (25) 



The next point is the computation of the extrema of 
the radial and polar motions. For the radial motion there 
are four extrema [see Eq. (90)]: r apo > r pcri > r 3 > r 4 . 
From Eqs. (19) and (90) we know that these extrema 



must satisfy the following relations: 

2M. 

r .po + r pc„ + r 3 + r 4 = YZTE2 ' 

'"apoVri + r 3 r 4 + ( r ap „ + V.-iX r 3 + r i) 
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2M.Q 

l-E 2 
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(E16) 
(E17) 
(E18) 
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From these relations we can find (r 3 ,r 4 ) using, for in- 
stance, Eqs. 1E151 and (|E18|), to get 



3. 4 = A ± VA 2 - B , 



(E19) 



where the plus sign corresponds to r 3 , and A and B are 
given by 



A = 
B = 



M. 

1-E 2 
a 2 C 



1-E 2 



(E20) 
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For the polar motion, and using the variable z = cos 2 9, 
there are two extre ma [ see Eq. (91)]: z + and z_, with 



z + > z_ 



From Eqs. (|20|) and (|9ip, these extrema satisfy 
E 2 ) 



+ z 4 



a 2 (l 



C 



a 2 (l-E 2 ) 
C 



a 2 {l-E 2 ) 



(E22) 
(E23) 



from Eq. (El), we can find z, from any of 



Given z_, 
these two. 

Let us now consider the inverse case in which we know 
the set (E, L Z ,C/Q) and we want to find the orbital pa- 
rameters (p, e, inc ) and other important quantities. We 
can start from the equations for the extrema for the radial 
and polar motions. In the first case, any extrema r+ will 



satisfy the following quart ic [from Eqs. (19) and (90)]: 

^ + 03^+03^ + 0^ + 00=0, (E24) 
where the coefficients are given by: 
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We can solve this quartic equation by steps (see, 
e.g. |115j V First, let us consider the following cubic equa- 
tion: 



y 3 + b 2 y 2 + b x y 



b = 0, 



where the coefficients are given by 



5 *. 
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1 „, 
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-Se 



26 2 
1 



(E27) 
(E28) 

2 2 8 ' (E29) 

and where S, r, and e are the coefficients of the depressed 
quartic u 4 + Su 2 +tu + e = . This equation is associated 
with the initial quartic via the change of variable: r + = 
u — a 3 /4. Then, the relations between the coefficients of 
the depressed quartic and those of the initial quartic are: 

3 n 1,1 



Now, let us consider a real solution to the cubic equation 
above, namely y 1 (a cubic equation with real coefficients 
always has at least one real root). Then, the four solu- 
tions of the initial quartic can be written as follows: 



\a 3 + g { s iV S + 2 Vi 



3(5 + 2y x + s 



2t 



1 V S + 2 Vi 



, (E32) 



where s x and s 2 are two independent signs, which lead 
to the four solutions. We can then immediately identify 
the extrema r apo > r pcri > r 3 > r 4 and from Eq. ( 22 ) we 
find e and p. 

In the case of the polar motion, things are a bit eas- 
ier. We need to find the roots of the following quadratic 
equation for z+ = cos 2 8 [from Eqs. (20) and (91)]: 



z\ + Cl ^ + c = 0, (E33) 

where the two coefficients are given by 

_ a 2 {l-E 2 ) + L 2 z + C 
Cl ~ "~ a 2 (l-E 2 ) 
C 

C ° ~ a 2 {l-E 2 ) ' 
Then, the two extrema of polar motion are given by 



(E34) 
(E35) 



2c n 



± 



2c n 



(E36) 



and from Eqs. (92) and (24) we find the inclination angle 



The inclination angle i follows from Eq. (25) as 



before. 



Appendix F: formulas for the Fundamental 
Frequencies and Periods 

In order to give formulas for the fundamental frequen- 
cies and periods with respect to the Boyer-Lindquist co- 
ordinate time, it is very convenient to start first consid- 
ering the frequencies and periods with respect to a new 
time (see [66]): d/dX = p 2 d/dr, which separates the ra- 
dial and polar dependence in the sense that it terms of A 
the equations for r and 9 become: 



R(r), 



d9Y 
dXJ 



6(0), (Fl) 



where R(r) and Q(9) denote the right-hand sides of 
Eqs. (fl9|) and (pOl respectively. The fundamental fre- 



quencies of the radial and polar motions associated with 
the A time are T r = 27r/A r and T g = 2ir/A e , where the 
periods, denoted by A r and A e , are given by 
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Then, following [57] , the radial and polar frequencies with 
respect to the A time can be written as 



T = 7iy(l - £ 2 )(r apo - r 3 )Q pcr . - r 4 ) 



where 



2K{k r ) 



iray/{l - E 2 )z^ 



( r , P o-Vi)( r 3- r 4) 

(r apo -7- 3 )(r pcri -r 4 ) ' 



(F4) 



(F5) 



and JC denotes the complete elliptic integral of the first 
kind (we adopt the definitions of |115j for the elliptic 
functions). 

The right-hand side of the equations for t and (f>, 
Eqs. (18 1 and (21), have the following structure in terms 



of the A time 



^L=T r {r)+T e {e) + aL zl 



dX 



(F6) 
(F7) 



and from here we can compute the azimuthal frequency 
with respect to the A time,r T^, and also the rate of 
change of the time t with respect to the time A, T t , as 
follows (see [87] , but note that there is a typo in the first 
expression in equation 21, our Eq. (F9), there is a closing 



square bracket that is located in the wrong place, the one 
that has E/2 as coefficient): 



2a T r 
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where £ and II denote the complete elliptic integrals of 
the second and third kinds respectively, and 



,)(r a 



(» 



:) 



(F10) 
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From these expressions for the fundamental frequencies 
and periods with respect to the A time we can compute 
the corresponding ones for the Boyer-Lindquist coordi- 
nate time t by using the following simple expressions: 



n, = 



T . 



T, 



(F12) 



and 



T = 



2n 



T, = 



2tt 



2tt 



(F13) 



Appendix G: formulas for the Evolution of Circular 
non-Equatorial Orbits 



Here we provide the expressions of the components of 
the matrix that determines the evolution of the Carter 
constant and radius of circular nonequatorial orbits in 
terms of the evolution of the energy and angular momen- 
tum component along the spin axis [see Eq. (36)]. The 
expressions of the c AB (A, B — 1,2) and d quantities, as 
functions of (A/., a; E, L z , C, r Q ) are: 
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c n = 4E(1 - E 2 )r 6 - YIM.Erl + 2E [a 2 (l - E 2 ) + 3(C + h\)\ r 4 a + 8M. [2a(L z - aE) - E(C + L 2 Z ) + a 2 E 3 ] r 3 a 
+ 2a [a 3 E(l - E 2 ) + aE(C + L 2 Z ) + 6M 2 (aE - L z )] r 2 Q - AM,a 2 E [(aE - L z ) 2 + C] r D 



+ 4aM 2 {L z - aE) [(aE - L z ) 2 + C] , (Gl) 

c 12 - 4(1 - E 2 )L z r 4 a - 16M.(1 - E 2 )(L Z - aE)r 3 Q + 2 [6M 2 (L Z - aE) -L Z (C + L 2 Z + a 2 (l - E 2 )} r 2 Q 

+ 4M.L Z [(aE - L z f + C]r Q ~ AM 2 (L Z - aE) [(aE - L z ) 2 + C] , (G2) 

c 21 = -2 [Erl - 3M.Ert + 2a 2 Er 3 a + aM,(L z - 2aE)r 2 Q 

+ a 4 Er Q - a 3 M,{L z - aE)] , (G3) 

c 22 = -2a [M.Erl - aL z r a - aM,(aE - L z )] , (G4) 

d = 2(1 - E 2 )r 4 - 4M.(1 - E 2 )r 3 Q + [6M 2 - C - L 2 Z + 5a 2 (l - E 2 )] r 2 + 2M. [a 2 (E 2 - 3) - 2aL z E + C + L 2 Z ] r a 

+ a 4 (l - E 2 ) + a 2 (C + L 2 ) - 2M 2 [{aE - L z ) 2 + C] . (G5) 
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